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ABSTRACT 
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CHi 



The calculation of non-lifting potential flow about arbitrary three 
dimensional bodies is examined in detail with specific interest in the XYZ 
Potential Flow program developed by the David W. Taylor Naval Ship 
Research and Development Center. The program uses a surface singularity 
distribution to solve the Neumann boundary value problem by means of a 
source panel method assuming a flat element with a constant source 
density over the area of the element. Boundary conditions are applied at 
control points on the elements producing a system of linear equations for 
the source density. When the source density is known, velocities and 
pressure coefficients may be calculated. 

The main purpose of this paper is to present the details of the 
approximation of an arbitrary three dimensional body using quadrilateral 
elements, and to provide a detailed derivation of the exact source panel 
integrations in order to gain insight for future research at Texas A&M 
University. A variation of the Hess method of surface discretization 
using quadrilateral source panels is described in detail as it is used in the 
XYZ Potential Flow program. The exact source panel integrations are 
derived in detail. 

A general discussion of other aspects of the program is included. 
Velocities and pressure coefficients for flow about a triaxial ellipsoid 
are calculated using the XYZ Potential Flow Program solution, and the 
results are compared with the analytical solution and the Hess Program 
solution. 
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1.0 INTRODUCTION 



This paper examines two aspects of the development of the XYZ 
Potential Flow Program (hereafter referred to as the XYZPF Program), a 
FORTRAN program which uses a source panel method to approximate 
solutions to steady potential flow problems about arbitrary three 
dimensional bodies. The aspects examined in detail are (1) the 
description of the details of the approximation of an arbitrary 
three-dimensional body using quadrilateral elements, and (2) a detailed 
derivation of the exact source panel integrations. 

The XYZPF Program was developed specifically for applications in 
numerical ship modelling and hydrodynamics studies at the David W. 
Taylor Naval Ship Research and Development Center (NSRDC) in Bethesda, 
Maryland. The format of the program is based on the work of Hess and 
Smith (1962) in the numerical calculation of non-lifting potential flow. 
A similar program is maintained by the Aerodynamics Division of the 
McDonnell-Douglas Corporation, referred to in this paper as the "Hess 
program.” The XYZPF Program is a modification of what has come to be 
known generally as the Hess Method. The most significant modifications 
are improvements to the method of solving the matrix equation for the 
source density, and greater flexibility in the input options (Dawson and 
Dean 1972). 

Though potential flow is a product of mathematics, and is never 
found in a real fluid, the results of potential flow calculations provide 
usable information for flow regions external to a thin boundary layer, 
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with little or no boundary layer separation. For such flow fields, the 
region outside the boundary layer may be considered to be effectively 
inviscid, and may be closely approximated by potential flow models. 

Small viscous effects can be accounted for by "thickening" the body by the 
appropriate displacement thickness. Displacement thickness accounts for 
the region of retarded fluid flow in the boundary layer inversely 
proportional to the square of the free stream velocity. Downstream of 
the point of boundary layer separation, the potential flow model no 
longer applies. 

Prior to the development of numerical methods, analytical 
solutions were generally restricted to simple analytical shapes (Kellogg 
1929). The need to solve boundary value problems for arbitrary 
boundaries in continuum mechanics has fostered the development of 
numerical approximations to the integral equation expressions. While 
the integration methods have been well known for quite sometime, only 
since the advent of high speed computers have many of the problems 
been practical to solve by numerical methods. Among the numerical 
methods being used are finite differences, finite elements, and the 
boundary element method. 

"Finite differences” and "finite elements" are numerical methods 
which satisfy the boundary conditions, and then approximate the 
solution to the governing equation in the fluid domain. These methods 
discretize the domain into a network of elements or cells. 
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Another approach is what is now known as the "Boundary Element 
Method," in which the governing equation is exactly satisfied in the 
domain, and the boundary conditions are applied through a boundary 
discretization method. The boundary value problem is reformulated as a 
boundary integral equation which is then discretized by subdividing the 
boundary into a finite number of surface elements. Each element is 
represented by an analytical function, and the source density function is 
integrated overthe surface of each element. Two factors governing the 
accuracy of the boundary element method are the boundary discretization 
method and the source panel integration. These two factors are examined 
in detail in this report, as a detailed derivation of the exact source panel 
integration, including the development of the source panel geometry, has 
not previously appeared in literature. 

The difference between the domain methods and the boundary 
methods is significant. The domain methods discretize the domain, while 
the boundary methods discretize the boundary. Thus, the boundary method 
reduces the dimension of the problem by one, as depicted in figure (1). In 
the application of the XYZPF Program, the problem is reduced from a 
three-dimensional problem in the domain to a two-dimensional problem 
on the boundaries. This method is well suited to problems in which the 
limits of the domain are infinite or difficult to define, in that the 
problem is applied to the boundary rather than the domain. 
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FINITE ELEMENT DISCRETIZATION 




BOUNDARY ELEMENT DISCRETIZATION 
Figure 1 . Discretization Methods 



Just as there are many variations of domain methods, there are 
also a variety of boundary methods. In general, they can be classified as 
"indirect” or "direct" formulations. The "indirect” method assumes a 
continuous source distribution over the surface of the body, and a 
solution which satisfies both the governing equation in the domain, and 
the boundary conditions on the body surface. The result is an integral 
equation on the boundary which has the surface source density function as 
its unknown. By enforcing the boundary conditions at control points on 
the surface, a system of equations is produced by which the source 
density may be determined. 

The "direct" method solves the velocity potential function through 
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an application of Green’s Second Identity requiring the solution of a 
source distribution and a dipole distribution on the boundary. The direct 
method has more physical significance to the general boundary value 
problem, and more versatility in its application as it can be applied to 
Neumann problems, Dirichlet problems, or mixed boundary value problems 
(Brebbia 1984). 

The simplicity and accuracyof the indirect method has made it 
attractive for many applications. The source panel method is an 
application of the indirect formulation of the boundary element method 
to the Neumann type of potential flow problem, for which the normal 
derivative of the potential function is prescribed on the boundary. 

1.1 OBJECTIVES 



The purpose of this paper is (1) to describe the details of the 
approximation of an arbitrary three-dimensional body using quadrilateral 
elements, and (2) to provide a detailed derivation of the exact source 
panel integrations for use in future investigations at Texas A&M 
University using panels of higher order geometries and source density 
functions. This paper is not intended to be a user's manual, though a 
general discussion of other aspects of the program is also included. 
NSRDC Report 3892 (Dawson and Dean 1972) is a summary of the XYZPF 
Program for those strictly interested in its use. 
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2.0 HISTORICAL DEVELOPMENT 



The foundations of the boundary element method were laid early in 
this century beginning with Fredholm in 1903 when he established the 
existence of solutions to the Neumann problem through a reconstruction 
of the problem using a discretized boundary (Kellogg 1929). The solution 
was determined to be the potential of a simple source distribution on a 
boundary with a continuous normal derivative for an infinite domain. 

Later works by Kellogg (1929) in potential theory demonstrated the 
application of the boundary integral equation method in electrostatics, 
heat transfer, flow in porous media, and fluid flow problems, but 
development was limited by the difficulty of obtaining analytical 
solutions. No significant advances were made until interest inboundary 
integral equation methods was revitalized with the advent of high speed 
electronic computers. Investigators were then able to discretize the 
boundaries and solve the integral equations numerically. This method of 
solution became known as the boundary element method. Early 
development of such numerical methods was pioneered by Hess and Smith 
(1962) and Jaswon and Symm (1963). Hess and Smith dealt primarily with 
the indirect formulation eventually leading to a solution for the three 
dimensional problem as described in this paper. In a parallel work, 
Jaswon and Symm developed a direct formulation approach to the two 
dimensional problem. The XYZPF Program is based primarily on the work 
of Hess and Smith. Hess has since developed a higher order panel method 
(Hess 1979) and Lefebvre modified the XYZPF Program for calculating 
velocity potentials for five degrees of freedom (Lefebvre 1982). 
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3.0 THEORETICAL DEVELOPMENT 



3.1 THE POTENTIAL FLOW PROBLEM IN THREE DIMENSIONS 



The governing equation for idea! (incompressible, inviscid, 
irrotational) flow is Laplace's equation: 



V 2 * = 0 



(0 



where $ is the velocity potential, and V 2 is the Laplacian operator. The 
XYZPF Program deals with steady, uniform flow of an ideal fluid about an 
arbitrary three dimensional body. The velocity components at any point 
within the flow field may be obtained from the negative gradient of the 
velocity potential, i. e. 



V 



-V* 



3$ 1 _ dfr j 

dx 3y 




( 2 ) 



The freestream flow is defined as a uniform stream of unit 
magnitude. 





V. 
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(3) 



The key to the boundary element method is the Divergence Theorem 
(Green’s Theorem) which relates a volume integral to an equivalent 
surface integral reducing the three-dimensional problem to a 
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two-dimensional one. The expression for Green's second identity is 
(Lamb 1924): 

J]J(*v 2 w - wv 2 $) da = JJ( > dr 

in which Q represents the integration over the three dimensional domain, 
andT represents integration over the two dimensional boundary. The 
partial derivatives are taken with respect to the outward normal, n. The 
weighting function, w, is usually chosen to be the fundamental solution 
for three dimensions, w = 1/(4Ttr), where r is the distance from the 
source to an arbitrary point on the boundary. 




Figure 2. Potential Flow in Three Dimensions 



Consider an arbitrary three-dimensional body with surface T, 
having an equation of the form F(x, y, z) = 0 where x, y, z are Cartesian 
coordinates of the global reference system as shown in Figure (2). The 
unit outward normal, n, at any point on the surface is given by the 
gradient of the function describing the surface divided by the magnitude 
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of the gradient, i.e. 



(5) 



±VF 

n = 

VF 



where the sign of the unit normal vector is chosen to ensure that the 
vector is an outward normal. The potential function $ describing the 
flow field must meet the following boundary conditions; 

a. V 2 $ = 0 (Laplace's Equation) (6) 

b. For an impermeable boundary, the velocity normal to the surface 
must be zero relative to the boundary (the Neumann boundary 
condition): 



' dn 'r 



0 



( 7 ) 



c. The velocity potential approaches the freestream velocity 
potential as the distance from the body goes to infinity: 



$ -+■ as r 



(8) 



The total potential at any point in the domain is composed of the 
freestream potential and the disturbance potential due to the body, 



$ = <Poo + 



(9) 
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The disturbance potential, 9, satisfies the following boundary 
conditions: 

a. V 2 <p = 0 (10) 

b. From equation (7), the velocity normal to the boundary due to 
the disturbance and due to the onset flow must be of equal 
magnitude, but opposite sign. Then from equation (9) 

(£) = n(p) • V„ (11) 

V On .'r 



Note that the normal vector is a function of position on the surface 
of the body. 

c. The disturbance potential approaches zero as the distance from 
the body goes to infinity, i. e. 

f -y 0 as | r | — ► w ( 1 2) 

3.2 MATHEMATICAL MODEL 



The disturbance potential of the body may be represented by a 
distribution of a source density function o over the body surface. The 
potential at an arbitrary point P(x, y, z) due to the surface potential is 
(Kellogg 1929): 
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(13) 



9 (x, y. 




cf(q) 

r(P,q) 



dr 



where q is the integration point on the surface, and 



r(P,q) = 



(x - x 0 ) 2 + (y - y 0 ) 2 + (2 - z*) 2 



is the distance from the field point P to the integration point q. 



The source density distribution function must satisfy the boundary 
conditions for the disturbance potential. Boundary conditions (10) and 
(12) are automatically satisfied by the form of the integrand. However, 
equation (11), the velocity normal to the boundary, combined with the 
Neumann boundary condition, equation (7), is the key to solving the 
boundary integral problem. 



The integrand becomes singular as the surface of the body is 
approached, i. e. | r | goes to zero. The singularity represents the local 
fluid flux normal to the boundary due to the local source density. The 
principal value of the singularity is -2tto(p), determined through a 
limiting process of the Gauss Flux Theorem (Kellogg 1929). The point p 
represents a field point which lies on the boundary. The integral 
expression is now composedof the contribution of the local source 
density and the contribution of the source density function over the 
remainder of the body surface. Solving for the velocity normal to the 
surface yields the following expression; 



1 1 



dr 
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From equation (II), this expression becomes: 



2tto(p) - fff to)' 
JJ 9n r(p,q) 



df = - n (p) 



00 
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This equation is a two dimensional Fredholm integral equation of the 
second kind, which ensures a unique solution, and that the diagonal 
elements of the system matrix will be dominant, each having a value of 
2 tt (Kellogg 1929). Once equation (15) has been solved for the source 
density a, the velocity components at any point of the flow field may be 
obtained by differentiating the disturbance potential function (13) with 
respect to the coordinate directions and adding the components of the 

freestream flow, V M . 



V(x, y, z) = V*. 



a<p i _ i 

ax dy J 




(16) 



The body shape does not have to be slender, axisymmetric, or 
simply connected, allowing for analysis of interior flow and a wide range 
of applications of the method. The only restriction imposed on the form 
of the body is that it must have a continuous normal vector. 
Discontinuities in the right hand side of equation (15) will produce 
unwanted singularities. Thus, in the process of approximating a body 
which has distinct corners, where there is clearly a discontinuity in the 
normal vector, the corner must be replaced by a surface with some finite 
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where is the static pressure at infinity. 



3.3 NUMERICAL MODEL 



In order to represent the surface of a body in the domain 
mathematically, the body may be described by analytical expressions 
which may provide an exact representation of the surface. However, the 
types of bodies which can be adequately described by such methods are 
severely limited. Another way to represent the body is to use a large 
number of analytical expressions, each describing only a small portion of 
the body. Hess and Smith (1962) suggested the use of an assembly of flat 
quadrilateral elements to model the actual surface of the body, as shown 
in Figure (3). Each quadrilateral approximates a region of the surface 
described by points which lie on the actual surface of the body. As planar 
elements, these quadrilaterals are clearly analytical, and when carefully 
constructed, the elements can approximate arbitrary three dimensional 
body surfaces without restriction. 




Figore 3. Approximation of the body by surface elements 
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curvature. However, application of this method has shown that the flow 
calculations give correct results for convex corners, while unrounded 
concavecornersmay or may not produce significant error, depending on 
the angle produced by the corner (Hess and Smith 1962). 

Because of the method of approximation, the calculation of flow 
velocities on the body surface are restricted to the points at which the 
boundary conditions were applied. Velocities at points other than those 
must be obtained by interpolation. Direct calculation of velocities at the 
edge of an element yields infinite velocities. 

With the solution of the system of linear equations for the source 
densities, the flow velocities at any point in the domain may be obtained 
from equation (16), and pressure coefficients are then computed from the 
velocities using a form of the Bernoulli equation: 



where P(t) is a constant independent of position. In the XYZPF Program, 
the flow is steady. Therefore, equation (17) can be reduced to 



and the pressure field can be expressed in terms of a dimensionless 
pressure coefficient C p as: 




p + -j- p | V = constant 



( 1 6) 




V 



2 



(19) 



V. 



00 
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The XYZPF Program uses the discretization procedure described by 
Hess and Smith (1962) with some minor modifications. The three 
dimensional body surface may be described using a large number of plane 
quadrilateral elements, each assumed to have a constant source density 
over the area of each element. Regions of the body requiring higher 
resolution for sharp curvature or anticipated velocity gradients will 
require a higher concentration of elements. 

Because the plane quadrilateral elements cannot fit edge to edge on 
a rounded surface, small gaps in the panel approximation contribute to the 
error of the approximation. However, the error due to the gaps is 
negligible when compared with the error of the basic model, that is, using 
flat panels to approximate a curved surface. Triangular elements have 
been suggested in an attempt to eliminate the gaps (Levy 1959), but the 
increased accuracy is so small that it may not justify the additional work 
of organizing the triangular elements in lieu of the simpler quadrilaterals 
(Hess and Smith 1966). The method presented is valid for an polygonal 
element with any number of sides. 

Equation (15) can now be decomposed into a summation of 
integrals, each representing the contribution of one element of the body 
surface. The unknown source density can be taken outside the integral, 
since it is assumed to be a constant over each element. The integration is 
performed over the area of the source element, and the boundary 
condition equation (11) is then enforced at a single point p in each 
remaining element. By performing this operation at each element of the 
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surface, a system of linear equations is generated which is equal in 
number to the number of surface elements and the number of unknown 
source densities. Equation (15) can now be approximated by the matrix 
equation (Dawson and Dean 1972): 



It is important to note that the influence coefficients Cjj andCjj are 

functions of geometry only, and once computed, need not be recomputed 
for analysis of several different flows. From the solution of equation 
(20) on the discretized surface, equation (13) may be applied at any point 
in the domain. Then, the velocity at an arbitrary field point P(x, y, z) in 
the domain may be determined from equation (16). With the velocity 
known, the pressure coefficient is determined from equation (19). 




(20) 



where 





16 



4.0 ORGANIZATION OF THE PROGRAM 



The XYZPF Program is actually composedof seven independent 
programs, referred to as sections PF1 through PF7. each of which builds 
on data generated from a previous section. This type of organization 
allows the user the flexibility of rerunning portions of the program using 
different flow parameters without having to go through the time 
consuming process of recalculating the influence coefficient matrix, 
which is dependent only on the geometry of the body. While the NSRDC 
program is very similar to the Hess program, there are also some 
significant differences. The following list of differences is taken from 
NSRDC Report 3892 (Dawson and Dean 1972): 

(1) The input to XYZ PF is arranged to facilitate the preparation of 
input for a series of problems in which only one part of the body is 
changed. Also, a number of checks are made on the input to help detect 
errors. 

(2) An option was added for the recomputation of the source 
density and velocities for only part of the body when only small changes 
are made. This option also provides for the use of the solution of one 
problem as an initial guess for the solution of another problem. 

(3) The matrix of influence coefficients is computed column by 
column instead of row by row. This column arrangement was used for the 
original LARC computer version because it required much less high speed 
memory. The computation is also about I0£ faster this way than with 
the row-by-row arrangement. 

(4) A simultaneous displacement iterative scheme is used to solve 
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the matrix coefficient for the source density. The scheme is slower than 
the successive displacement (Gauss-Seidel) scheme used in the Hess 
program, but it can be carried out using the matrix column by column 
instead of row by row. 

(5) When possible, an extrapolation procedure is used to reduce the 
number of iterations required for convergence. One such method is the 
Richardson extrapolation. 

The methods used in the XYZPF Program will be discussed in detail 
in the following sections. 
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5.0 DETAILS OF THE SURFACE APPROXIMATION 



5.1 PREPARATION OF THE INPUT 

Section PF1 is set up to read and store the input data, and to 
examine the cornerpoints of the quadrilaterals to detect obvious errors in 
the input. Because of the number of points which may have to be entered 
fora complex geometry, the user input is a major source of program 
error, and this first look for input errors will save lot of run time in the 
program as a whole. If Section PF1 detects major errors in the input, the 
program will not continue with the calculation of the coefficient matrix, 
but will stop and identify the grid location of the error. Minor errors may 
not cause the program to stop, but will be noted in the output. 

One of the major advantages to this program is in the organization 
of the input data. The surface is input in sections so that small portions 
of the input geometry may be changed without having to recalculated the 
points for the entire surface. The program also takes advantage of 
symmetry to minimize the input effort. Only the portion of the body 
which has no redundancy needs to be entered point by point. The 
remainder of the body is reflected across the planes of symmetry by the 
program to complete the surface representation. 

The surface is represented by a set of points in three-dimensional 
space which lie on the actual surface, and which will later be used to 
define the plane quadrilateral source elements. These points are defined 
in the global reference system. The points on the surface should be 
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selected in such a way as to provide an accurate representation of the 
surface with the fewest number of points possible. Portions of the 
surface which are highly curved will require a larger number of points to 
provide adequate resolution. Additionally, portions of the surface in 
which the flow field is expected to change rapidly will require a large 
number of points to accurately determine the flow field in that region. 
Some familiarity with fluid dynamics will provide a somewhat intuitive 
approach to properly distributing the elements. Elements should change 
gradually in size from areas of high concentration to those of just a few 
large elements, changing no more than 50 percent in size between 
adjacent elements (Hess and Smith 1966). The accuracy is only as good as 
that provided by the largest element in a particular area. The use of 
quadrilateral elements facilitates the use of known analytical equations 
and body contours to determine the input points. 




Figure 4. The 30 quadrilateral element in global coordinates 

For the purposes of this program, the body surface is represented 
by a large number of plane quadrilateral elements as shown in figure (3), 
each of which is assumed to have a constant value of source density over 
the area of the element. Each element is defined by four input points 
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which lie on the actual surface as shown in figure (4). Since each input 
point can be used as a corner for up to four elements, there is no need to 
enter the same point four separate times. The input points are organized 
in groups of four to form the quadrilateral element, and each point may 
also be associated with adjacent quadrilaterals. This is accomplished 
through the use of a two dimensional coordinate system in which the user 
assigns a pair of integers, m and n, to each point which identifies the 
"row" and the "column" in which it lies. A column of points will be given 
a common value of n, and each point in that column will have a unique 
value of m corresponding to the row in which the point lies. The 
orientation of these "coordinate" integers determines the direction of the 
outward normal for each element. Looking from the flow field toward 
the section of elements, if the values of m are increasing upward, the 
values of n must increase to the right. Increasing m and n can point in any 
direction with respect to the global reference system. In fact, the 
orientation can change from one section to another. However, any other 
relationship between m and n will produce an incorrect normal vector. 
Once assigned, the values of m and n also serve to identify the element 
for which the corresponding point is a corner. The four points which 
form a quadrilateral element are two points of one column, or n line, with 
consecutive m numbers, and two points of the next higher n line with the 
same m numbers as the previous two points. Thus, the element m, n is 
composedof the points identified by (m, n), (m+1, n), (m, n + l), and (m+1, 
n+1). 
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n=1 n=2 n=3 n=4 




Section 1 
CORRECT 



Section 2 
CORRECT 




Section 3 
INCORRECT 



Figure 5. Quadrilateral index numbers 



Each section of the body surface is formed by specifying a set of 
corner points corresponding to the m, n pairs for all of the quadrilaterals 
of the section. The user will sequentially assign an m number to the 
points for each n line, and also number the n lines for the section points 
entered. The first point in each n line will always have m = 1. The n lines 
are also numbered sequentially, but the value of n is not reset for each 
new section. The sequence of n numbers runs through all the sections as 
shown in figure (5). Points on a particular row or column do not have to 
be strictly colinear. By forming nearly triangular elements, a rounded 
planform can be approximated without conflicting with the numbering 
convention, as shown in figure (6). 
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Figure 6. Approximating a thin region with rounded planform (Hess & Smith 1 962) 



By entering data in sections, small changes in geometry can be 
performed without having to reenter all the points associated with the 
body surface. This feature is unique to the XYZPF Program, and offers a 
great deal of flexibility in design work. However, with the added 
flexibility comes more restrictions on both the input of the original 
geometry and on any modifications. There are four important 
restrictions on the input which are required to provide quadrilateral 
elements in groups of four to facilitate geometry calculations (Dawson 
and Dean 1972): 

(1) There must be an even number of elements in both the m and n 
directions in each section of the body. 

(2) The common corner point of a group of four elements must not 
coincide with any other corner point. The sides between the elements 
serve to define the local coordinate system, and serve as the axis of 
rotation when the surface is flattened for numerical differentiation of 
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the velocity potential. 



(3) Each set of four elements must have at least seven distinct 
corner points to allow the elements to more closely conform to a curved 
surface. This also allows for convergenceof N-lines orM-lines as may 
occur, for example, at the leading edge of an ellipsoid. Thus, only two of 
the four quadrilateral may degenerate into triangles by having two of 
their corner points coincide. This does not necessarily eliminate the 
possibility of more than two "triangular" elements since the adjacent 
sides of a quadrilateral may be colinear as shown in figure (6). 

(4) The normal vectors between two adjacent quadrilaterals in a 
group of four must be less than 90 degrees and preferably less than 45 
degrees. If a sharp edge is required, it should be a concavecornerwith 
respect to the flow field, and the input should be arranged so that the 
edge is along an outside boundary of the groups of four, and not through 
the center. 

When making small changes to the original geometry, the number of 
elements used in a new section must be the same as the number used in 
the original section unless the part being changed is at the end of the 
input data. Section configurations may be selected by natural divisions, 
as a matter of convenience to more easily handle large numbers of points, 
or as a tool to take advantage of symmetry. 

In setting up input data to use planes of symmetry, it is important 
to note that the XYZPF Program has certain restrictions on the choice of 
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symmetry planes. The user only has the option to select the number of 
symmetry planes. The planes which will be used as symmetry planes are 
preselected by the program to optimize the calculation procedure. 
Therefore, knowing this, the preparation of input data must consider the 
following restrictions imposed by the program (Dean and Dawson 1972). 

If only one plane of symmetry is used, the plane of symmetry is the 
y = 0 plane of the global coordinate system. As such, all the y 
coordinates of the input points must be of the same sign, i. e., all 
positive or all negative. If the body is closed and intersects the plane of 
symmetry, the points touching the plane, i. e., corresponding to y = 0, 
must also be entered with the input points. 

If two planes of symmetry are used, the planes of symmetry are the 
y = 0 plane and the z = 0 plane in the global coordinate system. Again, the 
y coordinates of all input points must have the same sign, positive or 
negative, and the z coordinates of all points must be of the same sign, 
positive or negative without regard to the sign of y. If the body surface 
intersects one or both of the planes of symmetry, the points which lie in 
the plane, i. e., those corresponding to y = 0 or z = 0, must also be entered 
with the input points. 

If three planes of symmetry are used, clearly the planes are the 
reference planes of the global coordinate system. As with the previous 
cases, all the x coordinates of the input points must be of the same sign, 
and similarly for the y and z coordinates. If any part of the body 
intersects any of the planes of symmetry, the points which lie in that 
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plane, i. e., x = 0, y = 0 or z = 0, must also be entered with the input 
points. 

5.2 SOURCE PANEL GEOMETRY 




Figure 7. The outer normal to the quadrilateral element 



With the surface points identified by the location numbers, m and n, 
and arranged in accordancewith program requirements, calculation of 
various aspects of the source panel geometry and formation of the plane 
quadrilateral element is the next step in the numerical integration 
process. Formation of all of the planar elements is identical, so the 
following discussion of source panel geometry will deal with only one 
characteristic element. The four corner points forming the basic 
quadrilateral are numbered in a clockwise direction from 1 to 4 as shown 
in figure (7). it does not matter which corner point is identified with the 
number 1 subscript, but the remaining points must be numbered 
consecutively in a clockwise direction when observed from the flow field 
in order to ensure an outward directed normal vector. These subscripts 
will be used to identify the cornerpoints for the remainder of this 
discussion. For this example, the points will be identified as follows: 
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Position Numbers 



Global Coordinates 



m+1, n+1 
m, n+1 



m+1, n 



m, n 




In forming the plane quadrilateral elements, the corner points, 
which are generally not coplanar, are used to form the local coordinate 
system, relative to the element. Recalling that the cross product of two 
vectors yields a vector solution which is perpendicular to both of the 
original vectors, the normal to the element may be obtained from the 
cross product of the diagonals of the element, 



where T 13 is the vector from point 1 to point 3, and T 2 « is the vector 
from corner point 2 to point 4. The unit normal is then: 



This unit normal now represents the first of the three local 
coordinate directions, this one in the C direction. The side of the 
quadrilateral from point 2 to point 3 is then used to obtain the second 
coordinate vector. 



N - T 2 4 x T 13 



( 21 ) 



n 




( 22 ) 
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n 




Figure 8. The second local coordinate vector 



N 2 - n x T 23 



( 23 ) 



and the unit vector 

n 2 




Similarly, the third local coordinate vector is obtained from the 
cross product of n 2 and n, the result of which is a unit vector. 



n,z 




Figure 9. The third local coordhate vector 



n 3 = n 2 * n ( 25 ) 

The unit vectors n 3 , n 2 , and n form an orthonormal basis and define 
the local coordinate system for the element in the T|, and C directions 
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respectively. Other methods of obtaining an orthonormal basis could be 
used just as well, and would make no difference to the remaining 
computations. The origin of the local coordinate system would most 
correctly be located at the "null point," the point at which the velocity 
potential has no contribution to the tangential velocity component on the 
source element. The null point is the point in each quadrilateral element 
where the normal velocity boundary condition is applied. However, with 
the exception of long, thin quadrilaterals, the physical difference between 
the null point and the centroid of the quadrilateral is not significant. The 
XYZPF Program will print a warning in the output when a quadrilateral is 
long and thin enough to jeopardize the accuracy of the approximation in 
that region. By locating the origin of the local coordinate system at the 
centroid, rather than at the null point, the difficult process of locating 
the null point for each element canbe eliminated, later calculations of 
the multipole expansion canbe simplified, and the boundary conditions 
canbe applied at the centroid without contributing significant error to 
the approximation (Hess and Smith 1966). Therefore, the origin for each 
local coordinate system is located at the centroid for the respective 
element. 
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5.3 LOCATING THE CENTROID 




Figure 10 . Locating the centroid of the quadrilateral 



The centroid of the element may be calculated by first dividing the 
area of the quadrilateral into two triangular areas, the triangles being 
separated by the line from point 2 to point 4. The area A, of the triangle 
defined by corner points 2, 3, and 4 is 



A, = 



2 



T 
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(26) 



Similarly, the area A 2 of the triangle defined by corner points 1, 2, and 4 
is 



A 2 = y |T,2 * Th | 



(27) 



In the global coordinate system, the X component of the centroid is 
given by 



where X] and X 2 are the averages of the X components of the corner 
points of each triangle. Substituting the values for X t and X 2 : 
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X 



■J At (X 2 + X 3 +X 4 ) + i'A 2 (X 1 + X 2 + x 4 ) 
A 1 + A 2 



1 

3 



( A 1 + A 2 ) X 2 + ( Ai + A 2 )X 4 + Ai X 3 + A 2 Xi 
Ai + A 2 
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x 2 + 



X 4 



Ai X 5 + A 2 X* 
A, + A 2 . 



Similarly 



Y 




Y d 



A1V3 + A 2 V t 
Ai + A 2 - 



Z = 



L-2 + 



A 1 Z3 + A 2 Zi 
A 1 + A 2 



(29) 



( 30 ) 



( 31 ) 



5.4 COORDINATE TRANSFORMATION 



Now that the local coordinate system is formed and properly 
located at the centroid of the element, the global coordinates of the 
corner points (X, Y, Z) are transformed to local coordinates (£, t\, 0 
through the components of the reference vectors of the local coordinate 
system as follows: 



f 




- ^ “ 




* 


n 3 x n 3y n 3z 




X - X 




£ 


n 2 x n 2y n 2z 




V - V 


- 


n 


n x Hy n z 




_ z - z _ 




vTV 

1 
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The corner points are projected into the plane of the quadrilateral 
element by setting the £ components to zero. The original diagonal 
vectors, T 13 and T 2 4 , will be on opposite sides of the resulting plane. The 
plane quadrilateral element is now completely defined. The program will 
sweep through all of the input elements using the assigned location 
numbers, and repeat this process for each element. 



5.5 MOMENTS OF INERTIA 

The calculation of the moments of inertia for each element are 
perf ormed f or use in the computation of the velocity coefficients using 
the quadrupole method. Any calculus text will give the moments of 
inertia of a planar section with a constant unit density about the origin 
to be: 




2 



Figure 1 1 . Forming the plane quadrilateral element 




(33) 



A 




( 34 ) 







A 
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For the triangular region defined by the corner points 2, 3, and 4, 
I XX = 77 (U 2 *i 3 ) 2 *U 3 *W 2 *U4^ 2 > J l (36) 

! a« = 77 I ( 1)2 + 113 > 2 + ( T13 + 3)4 ) 2 * < % ♦ % > 2 1 (37) 

A 

Ixy = -JJ U ^2 + $3 )( 'hs + '<13 ) + ( $3 + W ^3 + ^4 ) (38) 

+ ( £4 + ^2 )( TU + ^2 ) ) 



Similar equations can be generated for the triangular region defined by 
the corner points 1, 2, and 4. The moment of inertia for the entire 
quadrilateral is the sum of the corresponding expressions for each of the 
triangles. The resulting equations are: 



I XX = 77 IU 2 + ?3) 2 *U3*?4>M«4*«2> 2 1 

*7T l(«1 *i 2 } 2 *(«2 < ?4) 2 *(«4 + ?| ) 2 ) 



(39) 



A 



1 yy = J2 [ ( Ti2 + ^3 )2 + ( ^3 + T >,4 ) 2 + ( ^4 + ^2 > 2 1 < 40 > 

+ -jj i ( 'Hi + *n 2 )2 + ( ^2 + n '.4 ' ,2 + ( T l 4 + Tti ) 2 ] 



A 

Ixy = -JJ I ( S 2 + £?. )( T( 2 + T l3 ) + ( $ 3 + $ 4 )( T\l + ^4 ) 

+ ( $4 + ^2 )( ^4 + n2 )1 



+ “jy [( ^ + )( T l 1 + T| 2 ) + ( £>2 + ^4 )( + ^4 ) 



+ U 4 + *1 K T U + Tl t )] 



( 41 ) 
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6.0 THF MATRIX OF INFLUENCE COEFFICIENTS 



With the quadrilaterals completely formed, the next step is to 
calculate the velocities induced by the elements at the centroids of all 
the other elements. The total number of elements forming the surface 
will be represented by N. Let the source element be the (j)th element, and 
the element for which the velocity components are to be calculated at the 
centroid is the (i)th element. It does not matter how the (i)th elements 
are arranged in relation to each other as the sequence progresses. 
However, the sequence must be consistent as the calculations proceed 
from one source element to another. This program sweeps through the 
(i)th elements in the order of their location numbers, m and n. For each 
consecutive n line, the elements are swept in order of increasing m 
numbers. 

The result of the induced velocity calculations for the elements 
with unit source densities is an N by N square matrix of the values of 
induced velocities at each element, known also as the "matrix of 
influence coefficients." The XVZ potential flow program calculates the 
coefficients column by column, while the Hess program calculates them 
row by row. The advantage of one over the other depends on the method 
of later solving the matrix for the source densities. In calculating the 
influence coefficients, twenty-five quantities which describe the 
geometry of the source element are required to adequately calculate the 
induced velocity at the centroid of the (i)th element. These quantities 
include the coordinates of the centroid in the global coordinate system, 
the elements of the coordinate transformation matrix, the local 
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coordinates of the corner points, the maximum diagonal, the area, and the 
second moments of the quadrilateral element. Additionally, the Hess 
program uses the coordinates of the null point, making a total of 
twenty-eight quantities for that method (Hess and Smith 1962). 

When calculating row by row, the first (i)th element is selected, 
containing the "null" point, and the influence coefficients are computed 
for all of the (j)th elements in sequence before proceeding to the (i + 1)th 
element. This procedure requires the twenty-five quantities for each 
(j)th element to be available for calculation of the influence coefficient. 
Because each of the N (j)th elements is used N times with this procedure, 
calculating the geometric quantities or retrieving the values from low 
speed memory would be very time consuming, since the calculations or 
memory access would need to be performed N 2 times. Therefore, it is 
more practical to have the values available in high speed memory, 
although large matrices may exceed the storage capacity of high speed 
memory, imposing a limit on the number of elements which can be used. 
The advantage to the row-by-row calculation is that solution of the 
resulting matrix by the Gauss-Seidel reduction method does not require 
transposing the matrix, which would be another time consuming process 
(Hess and Smith 1962). 

Another alternative is calculation of the influence coefficients 
column by column. This method calculates the influence coefficients by 
sweeping all the (i)th elements for each (j)th element before proceeding 
to the (j+1)th element. Therefore, the twenty-five geometric quantities 
are retrieved from low speed storage only once for each (j)th element. 
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for a total of N times. This procedure is not limited by the capacity of 
high speed memory, and calculation of the coefficient matrix is 
approximately 10& faster than the row-by-row method (Dawson and Dean 
1972). This is the calculation method used by the XYZ Potential Flow 
Program. 



An influence coefficient represents the combined effects on one 
element of the velocity potentials of all the other elements comprising 
the body surface. For the quadrilateral element with a unit source density 
in the xy-plane, from equation (13), the potential at point P (x, y, z)due 
to the element is 





d£ dT\ 




+ (y - T|f + z 2 



n 



( 42 ) 



The integrand, 1/r, can be expanded in a series about the origin in 
powers of £ and T[. Each term of the series will contain the product of 
some powers of £ and T| with a corresponding derivative of 1/r 0 , where r 0 
is the distance of the field point P from the quadrilateral origin. 



and let 

1 



Then the series expansion through the second order term is (Hess and 
Smith 1962): 
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<P - AW - (M X W X + MyWy) + 1/2(1 XX W XX + 2 l X yW X y + lyy w yy) + • • • (43) 

The subscripts, x and y, used with w represent the respective partial 
derivatives. This series represents the multipole expansion of the 
velocity potential, since each term can be interpreted as a point 
singularity of a particular order. The first term is the potential at point 
P due to a point source of strength A located at the origin. The second 

term is the sum of two dipoles of strengths n x and n y located at the 

origin, oriented along the x and y axis respectively. The choice of the 
centroid of the quadrilateral as the origin of the local coordinate system 
causes the first moments, f1 x and My, to be zero. Therefore, the dipole 
terms dissappear, and are not dealt with anywhere in the program. The 
third term is the sum of three quadruples of strengths l xx , l xyi and l yy 

located at the origin. Kellogg (1929) shows that this secondorder 
approximation is absolutely and uniformly convergent, and Hess and Smith 
(1962) show that convergenceis rapid enough with an increase in r 0 that 
certain further approximations may be made without significant error at 
large distances r 0 from the source quadrilateral. 

Hess and Smith (1962) presented a comparisonof velocities 
calculated using the exact formulas, a simple point source, and a second 
order approximation. The comparisons were based on the ratio of the 
distance r 0 , between the centroid of the source quadrilateral and the field 
point P, to the length of the maximum dimension t, of the source 
quadrilateral, typically the maximum diagonal. The non-dimensional ratio 
is then r 0 /t. The results show that sufficient accuracy is maintained 
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while using a simple point source at ratios of (r 0 /t) > 4, using the second 
order source and quadrupole solution for the range 2.45 < (r 0 /t) < 4, and 
using the exact solution for ratios of (r 0 /t) < 2.45. In any case, the error 
goes to infinity as the field point approaches the edge of the quadrilateral 
where calculations indicate an infinite velocity. The XYZ Potential Flow 
Program uses a monopole source for (r 0 /t) > 4, the source - quadrupole 
formulae for 2 < (r 0 /t) < 4, and the exact formulae for (r 0 /t) < 2. Hess 
and Smith (1962) reported a maximum error of 0.001 in approximating any 
velocity component using the above criteria. 
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7.0 OFRIVATION OF THE EXACT SOURCE PANEL INTEGRATION 



From equations (2) and (42). the components of the velocity at the 
field point P(x. y. z)due to the source quadrilateral are: 



v v 




v y 



0f 

9y 



> 

Z 





(x - Q dg d q 
r 3 



(y - T[) d| dm 

F 3 




(44) 



(45) 

(46) 



Equations (44), (45) and (46) are evaluated by expressing each of 
the integrals as the sum of four terms, each term representing the effect 
of one side of the quadrilateral (Hess and Smith 1962). This method can 
also be generalized for polygonal elements with any number of sides. The 
potential function for each side of the quadrilateral is the combined 
potentials of semi-infinite strips whose boundaries are the side of the 
quadrilateral and two semi-infinite lines parallel to either the x or y 
axis. When observed from the domain, and the sides are traversed in a 
clockwise direction, the source strip on the right will have a source 
density of O = +1/2 and the source strip on the left will have a source 
density of o = -1/2 as shown in figure (12). When the sides are 
recombined to form the quadrilateral, the source densities outside the 
quadrilateral cancel each other, and the source densities within the 
quadrilateral combine to form a source density of o = +1. This will be 
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true for a planar element with any number of sides and in any relative 
orientation within the plane. 




7.1 THE Y VELOCITY COMPONENT 



From equation (45), the velocity component Vy is found by summing 

the four terms representing the contributions of the sides of the 
quadrilateral. For the side from point (£ t , T},) to point (£ 2 . t[ 2 ), the 
contribution is expressed as the integral over the area of the 
semi-infinite strips with the source densities of O = +1/2 and O = -1/2 
rather than the unit source density of equation (45). 



V u 






V = | f ' 52 « f f" 2 - f" ] 

J LJ-co Jti |2 J[ ( k _ £)* + {y - ^) e + 2 £ 
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Integrating with respect to r\: 



v « 12 = 2 { 



•7 2 ik' 



•[ (x - £) 2 + (y - T[) 2 + z 2 j 



[ (x - £) 2 + (y - T() 2 + z 2 fi 



Ttl, 



1 



-00 [ (x " £)* + (y - T|) S + Z 2 ]' ^ 



DO 



1 






"^12 



} 



The terms evaluated at t\ = +<*> and Tj = -oo cancel, and the terms evaluated 
at t\ - H 12 add to obtain the following expression; 



V, 



i r* 



2 61 



V 12 



2 i [ (x - £) 2 + (y - ti) 2 + z 2 3 
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v «. a = jr - 7 - (48) 

J -A 

Equation (48) is changed to a function of arclength s by the relation 



d£ 



?•_ - ^ 
s 1 



*1 



\| U 2 -Si ) 2 + ( t i 2 - V 



(49) 



where d 12 is the length of the side of the quadrilateral from (£,, -q,) to 
(£ 2 . ^ 2 ) as shown in figure (13). 
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^ * 


H i n 



Figure 1 3 . The potential due to a finite line source (Hess & Smith 1 962 ) 



Substituting equation (49) into equation (48) 



From figure (13), it can be seen that, in terms of arclength s, the distance 
r from point P to any point on the line from point 1 to point 2 is given by 



r = 4 f 2 + (si - $)- 



Substituting into equation (50) yields 



v y 1£ = 



i^ - /- 

s2 - cj r 

^12 Jo 



ds 



N 



f 2 + (Sj - s ) 2 



K2 - Ki f ds 

^12 J 0 4 f 2 + (S - Si) 2 



(51) 
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Evaluating the integral 



V, 



*2 -Si 



yi£ - 



log 



(s -sO + a] f2 + (s - St } 2 



d l£ 



S2 - Si 



^{iog[(d 12 -s,)*j f2t(dl2 . S]) 2 



- loq 



(-SiKjf 2 + (-s,) 2 



$2 - £1 



-^7— [log(r 2 - s 2 ) “log(r, - s t ) \ 



S2-S1 , (r 2 -s 2 ) 

T^ log T^y 



(52) 



The quantities r lt r 2 , s 1( and s 2 used in equation (52) are as shown in 
figure (13). Equation (52) is singular when r, = s 1t which occurswhen the 
field point P is located anywhere along the line defined by the side of the 
quadrilateral. This singularity may be removed by using the law of 
cosines (Hess and Smith 1962). 
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From equation (52) 



V, 



y i£ - 



2 — £1 , r 2 (i - cos £ 2 ) 



log 



(53) 



di2 n ( 1 ~ cos ) 

where 3, and fJ 2 are the interior angles shown in figure (13). Applying the 
law of cosines to figure (13) 



COS $| 



If + d| 2 2 - r 2 2 
2 r, dis 



(54) 



COS $2 



r 2 2 + d | 2 2 -rf 
2r 2 d| 2 



if - di 2 2 - ri 

2r 2 di2 



(55) 



From equations (54) and (55): 



r 2 (i - cos &>) 

HO - cos £i ) 



To 



rf - di 2 " - r 2 " 
2 r 2 d |2 



if + d, 2 ‘ 



r, 2 



2 r.d| 2 



Or 



r 2 d| 2 - rf + d| 2 c + r 2 
2r 2 d|2 






2 r.di 2 - r, 2 - di 2 * + r 2 £ 



2 r, d| 2 



(r 2 + d l2 ) 2 - f 



2r 2 d| 2 - f + d )2 2 + r 2 2 _ 

2 r 1 d, 2 -r t 2 -d, 2 2 + r 2 2 " -(r, - d, 2 ) 2 + r 2 2 



1(1*2 + d| 2 ) + r |l 1(1*2 + d| 2 ) “ f *| ] l’| + 1*2 + d| 2 

( r 2 + (r, - d| 2 ) J (r 2 -(r, - d J2 ) ] r i + r 2 - d i2 



(56) 
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Substituting equation (56) into equation (53) yields the final form of the 
exact equation of the y component of the velocity induced by the side of 
the quadrilateral from point 1 to point 2: 




£2 — £l r | +r 2+d|2 

— -r loq 

d 12 r i + r 2 - d i2 



(57) 



Equation (57) is applied to the remaining sides of the quadrilateral 
simply by substituting the appropriate point numbers for the corner 
points of each side. The total contribution of the quadrilateral to the y 
component of the velocity is the sum of the four terms representing the 
contributions of each of the sides. The y component of the velocity at 
the field point P is now given by: 



V u = 



■ s 2 - C] 



'12 



, r ,+ r 2 + di 2 $3- $2 , r 2 + r 3 + d„ 

l0 9 T— 7 — + — 1°9 

sJ 



r, + r 2 - d| 2 



r 2 + r 3 - d 23 



+ 



C4- e 3 , r s + r 4 + d, 4 
— q loq — 

a 34 f 3 + f 4 " d 34 




log 



r4 + n + d 41 



r 4 + n - d 4J 



(5S) 



7.2 THE X VELOCITY COMPONENT 



A similar derivation process is used to produce the equation for the x 
component of the velocity induced by the side of the quadrilateral from 
point 1 to point 2. The semi-infinite source strips are constructed 
parallel to the x axis, and the order of integration is reversed. 
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The x component of the velocity at the field point P due to the 
quadrilateral is given by: 



ite-Tli 


- log 


r l + r 2 + d | 2 1 f l 3 - T l 2 


log 


r 2 + r 3 + d 2 3 


d 12 


r, + r 2 - dj 2 0 2 3 


r 2 + r 3 - d 23 


-' 0.3 


log ■ 


^ 3+^*4 + d 34 'Hi — ^4 


log 


r 4 + r, + d 41 


*34 


r 3 +r 4 -d 34 ' d41 


f 4 + n - d 41 



(59) 



7.3 THE Z VELOCITY COMPONENT 



The z component of the velocity at the field point P due to the 
quadrilateral is obtained in a similar fashion, using semi-infinite source 
strips, this time parallel to the y axis. From equation (46), the 
fundamental velocity potential of the semi infinite sourcestrips is 
integrated in a manner similar to that used to obtain equation (47), and 
the z component of the velocity due to the side from(£,, ti,) to (£ 2 , ti 2 ) is 
given by 



V, 
z 12 



« [i 


P - 1 


T| L - 




c r i 


r T u 2 r 
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zd'n 
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( 60 ) 



Performing the integration with respect to T|, the integral 

dm 



I 



[ (x - c) 2 + (y - m ) 2 + z'-'P 
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fits the integral form 



-J 



dF 



[C 2 + F 2 jn = 2C 2 (n-1) [(C 2 + F 2 ] n_1 



+ ( 



2n-3)f 

J [C 2 + I 



F 2 ] 



n-1 



where 



C 2 = (x - O 2 + z 2 



F = (y - ti) 

dF = -dT[ 
n = 3/2 



Then, from equation (60) 






-if 



S .2 



ai: 



, f 



(y - ti) 



c ,i 



'' l[(x - £;) 2 + z 2 H (x - $) £ + (y - + z 2 j 

(y - 'a) 

[(x-C ) 2 + z 2 ]I(x-$) t + (y-n ) 8 + z*] )s 

(y - j 0 

Kx - O 2 + z 2 ][ (x - O 2 * (y - r\ f * z 2 ]^ 

(y - ti) 

[0<:-£j 2 + z 2 3[(x-a 2 + (y-T l ) 2 + z 2 ]^ 



O 1 



T l|£ 



-CO 



co 



(6 ! ) 



*n.|2 ** 



Again, the terms evaluated at +oo and -oo cancel and the terms evaluated at 
ri 12 add to obtain the following expression: 



v =, 2 = - z { 



^2 



(y " Tlia)d£ 



s, [(x - £) 2 + z 2 ][ (x - £) 2 + (y - 



'q, 2 ) 2 + z 2 



(62) 
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Without a convenient substitution with which to integrate equation (62), 
the integration is performed directly. Recognizing that along the line 
defined by the side of the quadrilateral from (£,, tj,) to (£ 2 , *n 2 ), t \\ 2 may 
be expressed as a function of £: 

= m 12 £ + b 12 (63) 

where the slope of the side, m 12 is given by 

~ (64) 

hii2 = ^ * 

£2 - £i 

and b|2 may be determined knowing that T| 12 = t \\ when £ = 

6 T li-^'n£ (65) 

bl2 = 

Substituting equation (63) into equation (62) yields 



„ _ _ f (y -m, a ^-b| g )d^ 

vz 12 " J i, l(x - O 2 + z £ ][ (x - O 2 * (y -m !2 C - b l2 ) 2 + z 2 {b } 



Define the quantities 



^12 = y “ b 12 " m 12 x 
u = x - £ 



(67) 

(68) 



- Then du = -d£ 



y - b 12 - m 12 fc = m 12 u + q 12 



(69) 

(70) 
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By a change of variable, equation (66) is now expressed as a function of u: 






V; 



. -. z f 



(m I2 u + q, 2 ) du 



J „ r [ u 2 + 2 2 ] [ u 2 + ( m 12 u + q, 2 ) 2 + z 2 



(71) 



, (m| 2 u*g, g )du 

J [u 2 + z 2 ][(m I2 2 + i)u 2 + 2m l£ q l2 u + q l |+ z 2 ] 

X- ci 



which fits the form of 



I 



(Lu + M) du 



(Au 2+ 2Bu + C)4 (au 2 + 2bu + c) ' 



(73) 



where 



L = m 12 


n = q, 2 


A = 1 


a = (m 12 2 + 1) 


o 

II 

CD 


b - m 12 q 12 


C = z 2 


c = q 12 2 - z 2 



From Hardy (1944), this integral form may be integrated by the 
substitution 



where ji and u satisfy 

ajiu + b(p + u) + c = 0 
Ajiu + B(ji + u) + C = 0 
- and are the roots of the equation 

(aB - bA) £ 2 - (cA- aC) $ + (bC - cB) = 0 (77) 



(75) 

(76) 
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Substituting the appropriate values into equation (75), the roots of the 
quadratic equation are 



J-i = - 



Jia. 

m ( . 



(78) 



v 



mi 2 z 2 

^12 



(79) 



It can be verified that these values satisfy equations (75) and (76). 



Substituting equations (78) and (79) into equation (74) 



m i £ 2 2 q 12 t 

qia " m is ( 80 ) 

t + i 



du = 



qiat 

mi 2 



rn 1^2' 

“v” 



(t + I ) 2 



dt 



(8 1 ) 



By substitution and a change of variable, equation (72) becomes a function 
of the parameter t. After simplification, the integral now fits the form 
of 

r dt 

K J (c<t 2 + £)4(tft 2 + 5)' {62) 



where 



K = -(m 12 6 q 12 z 4 ♦ 2m ]2 4 q, 2 3 z 2 + q 12 5 m 12 2 ) 
o< = q 12 4 ♦ m 12 2 q 12 2 z 2 



50 



3 = m 12 4 z 4 + m ]2 2 q 12 2 z 2 

# = Qi2 4 + m, 2 2 qi2 2 2 2 

8 = m 12 6 z 4 + m 12 4 z 4 + 2m, 2 4 q 12 2 z 2 + m 12 2 q 12 4 + m 12 2 q J2 2 z 2 
Equation (82) can be rationalized by the substitution 



v = 



t 



s 



from which it can be shown that 



(83) 



t £ = 



v 2 & 



1 - V 2 '6' 



(64) 



dt = 



(I - vW 



% 



dv 



(85) 



Substituting equations (84) and (85) into equation (82) and simplifying 
yields the integral in terms of the parameter v: 

, r dt r dv 

k J (<x t 2 + + 6)' " K J 3 + (« S - 3 3 ) v 2 0 ■' 



which fits the form 

dv 

a 2 + b s v* 



J 



1 . -i bv 
tar > — 



where a 2 = 3 

b 2 = (cxS - £20 



(87) 
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Performing the integration 



K 



r dv 


K 








j $ + (oc 6 - $ Z 


>)v £ 


■y lan 

j 


V 

a 


£ J 



(88) 



From equations (80), (83), and the expressions for «, £, 5, and from 
equation (82), and after a considerable amount of algebraic manipulation 
and simplification, equation (88) becomes 




tan 



y 


c< 6 “ ■£> o 


•V 


— i 



4 - tan -1 



m 12 z 2 - q 12 u 



z Jz 2 + u 2 + (m 12 u + qia) 2 ' 



( 89 ) 



From equations (63) and (67), equation (89) becomes 



tan 






m 12 z 2 - Q 12 u 



2 z 2 + u + (rn, £ u + q 12 ) 2 



1 

— tan 



m 1 £ (u 2 + z 2 ) - (y - T| 12 ) u 

z J u 2 + (y - T( 12 ) 2 + z 2 ' 



( 90 ) 
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Finally, applying these results to equation (71) 



V-T 

x 



r tm| 2 u + q l£ ) du 

Z|2 ~ J . [u 2 + z £ ][u 2 + (m, 2 u + q i2 ) £ + z 2 ]— 

X Sj 

rn 12 (u 2 + z £ ) - (y - r^) u 



= -tan 



z Ju 2 + (y - T[ 1£ ) 2 + z 2 ' 



x-z,. 



**S| 



- tan' 



-1 



m 12 ((x - ) a + z 2 ) - (y - t\ 12 )(x - £i ) 
z -J(x - Ci ) 2 + (y - T[ 12 ) 2 + z 2 ' 



tan' 1 



m 12 ( (x - & ) 2 + z 2 ) - (y - ti 1£ )(x - 6 ) 
z 4(x - 6) 2+ (y - iq 12 ) 2 + z 2 ' 



(91) 



Recall that when x = £ 1( y = *q 1 and when x = $ 2 , y = ti 2 . Then, for the sake 
of a more compact equation, define the following quantities: 
e, = (x - W ♦ z 2 e 2 = (x - £ 2 ) 2 + z 2 

hi = (y - t}i)(x - $,) h 2 = (y - t| 2 )(x - $ 2 ) 



The quantities r, and r 2 are as shown in figure (12), where 



n = 4(x - $i) 2 + (y - rt,) 2 + z 2 ‘ r 2 = J(x - $ 2 ) 2 + (y - t| 2 ) 2 + z 2 ‘ 



Substituting these quantities into equation (91) yields the exact z 
component of velocity due to the side from point (^/q,) to (£ 2 /q 2 ) in the 
form used by the XYZ Potential Flow Program: 



, r m 12 ei - hi 

V 7 = tan 

^12 |_ Z 



-tan' 1 



m 12 e 2 - hg 
z r 2 



(92) 
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The total z component of the velocity at the field point P(x, y, z)due to 
the quadrilateral element is the sum of the four sides: 



= tan' 1 


FT! i o 6 1 

^ n 


- hi 


• 

i 


rni2 ~ hs 
. 2 


+ tan' 1 


ffi o 2 


- h 2 ' 


-tan' 1 


rfioj “ h 


z r? 




. 2 r 3 


+ tan' 1 


' m 34 e 3 
2 r s 


*“ hj 


-tan' 1 


^3464 ■ h4 

. 2 r 4 


+ tan' 1 


m 41 t: 4 


- h 4 ' 


-tan' 1 


rn 41 e 1 - hi " 


. 2 r 4 




2 ri 
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8.0 APPROXIMATIONS OF THE INDUCED VELOCITY 



8.1 QUADRUPOLE METHOD 

As previously mentioned, as the ratio of (r 0 /t) exceeds the value of 
2, then certain approximations may be made which greatly reduce the 
calculation effort otherwise required by the exact method. In the range 
of 2 < (r 0 /t) < 4, the XYZ Potential Flow Program uses the second order 
approximation of the potential described by equation (43). With the 
origin at the centroid of the quadrilateral, the first moments are zero, 
and the second order approximation is 

<P = AW * (l/2)(l xx w xx * 2! xy w xy * I yy w yy ) (94) 



where the first term is a point source of strength A, the second term is 
composedof three quadrupoles of strengths l xx . l X y. and lyy located at 

the local origin, and the subscripts on w indicate the partial derivatives 
of w with respect to those variables as before. The quantity A is the area 

of the element, and the terms l xx , l X y, and lyy are the respective moments 

of inertia of the source element given by equations (39), (40), and (41). 

To obtain the velocity components at the field point, equation (94) is 
- differentiated with respect to the coordinate directions giving: 
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(95) 



V x = ~ [ Aw x + 2 ^xx W xxx + IxyW XX y + 2 I y y W X y y J 

Vy = - = -[AW U + \ Ixx w xxy + IxyWxyy + ^IyyWyyy] (96) 

V z = - = - [ Aw z + \ I xx w xxz + IxqW xyz + \ IyuW yyz ] (g?) 

Recall inq that w = . I — = -J- 

y J x 2 + y 2 + z 2 r o 



the derivatives of w, as expressed by Hess and Smith (1962) and as used in 
the XYZPF program, are 



Ww = 


-x r 0 -3 








w y = 


— p r “3 

y r o 


\ 

s 




(99) 


w z = 


-zr 0 -3 ^ 








w 

XXX 


= 3x(3p + 


10 x 2 )r 0 -7 






W X xy 


= 3y p r c -' 








w 

‘xyy 


= 3xqr 0 -' 








Wy y y 


= 3y(3q + 


1 0y 2 ) r 0 -7 


> 


(100) 


w xxz 


= 3z p r 0 w 






W X ys 


= -15xyz 


w 1 






Wyyz 


= 3z q r 0 -‘ 


r 







where 



p = y 2 + z 2 - 4x 2 
q = x 2 + z 2 - 4y 2 
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8.2 MONOPOLE METHOD 



When the ratio of (r 0 /t) is greater than 4, then the quadrilateral 
may be approximated by a simple source corresponding to the first term 
of equation (43). Then the velocity components at the field point due to 
the quadrilateral are given by 



V v = -4 






= - HYYv 



tlx 

v »= =- Aw ' 

v - a,f ' 
v = - "IF AWj 



( 101 ) 



( 1 02) 



(103) 



where the partial derivatives of w are those given in equation (99). 
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9.0 SOI V1NG THF MATRIX EQUATION FOR SOURCE DENSITY 



9.1 JACOBI'S ITERATIVE METHOD 



From equation (20), the matrix equation may be solved for the 
constant source density o j for each element which satisfies the boundary 

condition equation (11). Equation (20) suggests the use of Jacobi's 
iterative method of matrix solution in the form 



(ni+1) 

Si 



Vi + £ Cij Oj (m) , 1=1,2 N 

j = t 
i* i 



(104) 



where N is the number of elements composing the body surface, and m is 
the number of iterations completed. A partial sum of equation (104) is 
computed for each of the ith elements before proceeding to the next jth 
element. The iteration is complete when the summation of equation 
(104) includes all of the jth elements. Because the values of the source 
densities at all of the elements are recomputed before any of them are 
used in the iteration, this method is also called the simultaneous 
displacement method (Ralston 1965). This is contrasted with the 
Gauss-Seidel iterative method used in the Douglas program. In the 

Gauss-Seidel method, as each new a j is computed, it is used immediately 

in the iteration process for calculation of . This is also known as 
the successive displacement method and is expressed as 
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(105) 



(m+1) 

<*i 



Vi ♦ 



1-1 

z 

j=1 



r w (m+1) 
Cij cJj + 



N 

z 

j=i+l 



CijOj 



(m) 



i = 1, 2, ... , N 



Though the Gauss-Seidel iterative method is faster, the Jacobi iteration 
method was selected for use in the XYZPF program in order to be able to 
perform the iterations column by column, since the coefficient matrix is 
also computed column by column, and the matrix does not have to be 
transposed for solution. 

When the (m+1)th iteration is complete, the values of the source 
densities are compared with those of the (m)th iteration and the 
differences summed for all of the elements. The total difference between 
successive iterations is then compared to a convergencecriteria input by 
the user. If the difference is less than the convergencecriteria, then the 
matrix solution is complete and the values of the source densities are 
stored for later use in computing velocities and pressure coefficients. If 
the convergencecriteria is not met, then the iteration process is 
repeated. After every five iterations, if the convergencecriteria is still 
not met, then an extrapolation is attempted in order to accelerate the 
convergence. The XYZ Potential Flow Program uses a Richardson 
extrapolation method, a numerical procedure which uses two approximate 
results to obtain a third approximation which is closer to the exact 
solution (Ralston 1965). 
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9.2 RICHARDSON EXTRAPOLATION 



The Richardson extrapolation assumes that the iterative process is 
convergent. For the iterative solutions S 0 , S, ( and S 2 , where S 0 is the 
most recent approximation and S 2 the oldest, the solution is convergent 
if 



While a Richardson-type extrapolation can take may forms, the XYZPF 
program uses a procedure developed from the following approximations 
(Dawson and Dean 1972). If there is only one dominant eigenvalue and a 
sufficient number of iterations have been completed, the iterative 
solutions may be approximated by 



Define the linear combination which, from equation (107), may be 




( 106 ) 



S 0 * S f + E X m 
Sj ~ Sj- + E X m_1 
S 2 ~ Sf + E X m " 2 
Sj * S f ♦ E X m-i 



(107) 



where S f is the true solution 



X is the eigenvalue 

E is the eigenfunction 

m is the number of completed iterations 
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approximated as 



A S 0 + 


(1 - A) Sj * S f + E X n_1 (AX + 1 - A) 


(108) 


The value of 


A may be chosen such that 




(AX ♦ 


1 - A)= 0 


(109) 


Then, from equations (108) and (109) 




AS 0 * 


(I -A)S,*S f 


(110) 



where the expression on the left convergesto the exact solution. 
From equations (106) and (109) 



Solving for A, 



Sg - 2 Si + S 2 0 

Since the value of A generally changes from element to element, a 
weighted average of A is used in the extrapolation, where 

N 

£ (S 2 (i)-Si(i)) (sign of D(i)) 

A= ^ ^ 

£ D(i) 

1=1 

Equation (113) is recomputed after every fifth iteration. If the difference 
between the new value and the old value is less than 0.02, then the 
solution is extrapolated. From equation (110), the extrapolated solution 
is 

S* = A S 0 + (1 - A)S, (114) 
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When there are two dominant eigenvalues, then the iterative solutions 
may be approximated by 

S( * Sf * E, X, m -' * E 2 X/"'' 

where Sf is the true solution 

Xj and X 2 are the eigenvalues 

Ej and E 2 are the eigenfunctions 

m is the number of completed iterations 

Define the linear combination which, from equation (115), may be 
approximated as 



B 2 Sq + B] + (1 ~ Bj ■ B 2 ) S 2 

~ S f + E, X/"" 2 [B 2 X, 2 + B, X 1 + (1 - B, - B 2 )] (116) 

+ E 2 X^ 2 [B 2 X 2 2 + B^ X 2 + (l ~ Bj ~ B 2 )] 

The values of B, and B 2 may be determined for which the eigenvalues Xi 
and X 2 are roots of the quadratic equation 

B 2 X 2 +B,X+(1 -B,-B 2 ) = 0 (117) 

Then, from equation (116) 

B 2 5 0 + B, S, + (1 - B, - B 2 ) S 2 * S f (118) 
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where the left side of the equation (118) convergesto the exact solution. 
Using equation (115) and eliminating terms containing E 2 : 



(Sq ~ S]) - X 2 (S] ~ S 2 ) - E] X^^Xi X 2 )(Xj - l) 

(5, - S 2 ) - X 2 (S 2 - S 3 ) = E, ~^2)(h ~ 0 (119) 

(S 2 -S 3 ) - X 2 (S 3 - S 4 ) = E, X 1 rr4 (X 1 -X 2 )(X, - 1) 

Solving for Xj 

^ 0 Si _i Xj i t Oj ” ) ( '— * 1 ~ } X c! V 2 “ O 3 ) ( ] ^ 0 

= (Si - Ss ) ~ Xs (Si- - S 3 ) = (S 2 ~ S 3 ) ~ X 2 (S 3 ~ S 4 ) 

From equations (117) and (120) 



(S4- S3 ) (So - 2 Ss +54) -(S4 - S2 ) l(Si - S2MS3- S 4 )] 

Bi= ( 1 2 1 j 



(S4- Ss ) (S4 -2S3 +32 )-(S 4 - Ss)[(Sl - S2 )— (S3 — S4)] , , 

Bs = ( 1 21) 



where D = (S 4 - 2S 3 - S 2 )(S 0 - 2S 2 + S 4 ) - (Sj - S 2 - S 3 + S 4 ) 2 



The weighted averages of B, and B 2 are used for the extrapolation as done 
with Ain equation (113). If the sum of the absolute values of the 
weighted averages of Bj and B 2 changes by less than 2%, then the 
extrapolation is performed. Then from equation (118), the extrapolated 
solution is 



= Bs So + B 1 S 1 + (1 -B, -B2)52 



(123) 
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10.0 CALCULATION OF VELOCITIES AND PRESSURE COEFFICIENTS 



With the influence coefficients and the source densities 
determined, the calculation of velocities is a relatively simple matter. 
From equation (9), the total velocity is the sum of the freestream 
velocity and the disturbance velocity due to the body. The product of the 
source densities and the influence coefficients are summed for all of the 
elements, and then added to the freestream velocity to determine the 
total velocity at any point in the domain. Velocities on the surface of the 
body are calculated at the null points only, as the boundary conditions are 
enforced only at the null point of each element, and velocities at other 
points in the element would produce significant error due to the method 
of approximation. The components of the velocity at the centroid of the 
ith element are 



Vi* = V, 



x 



V* - V- - 



; 1 - 



- v 

'’ 00 , 



* Z Cij x Cj 

j=1 

N 

+ Z Cij g Oj 
1=1 

N 

* Z Ci.i : CT) 

j=1 



(124) 



From equation (15), the velocity induced by an element at its own null 
- point has a magnitude of 2tt directed along the outward normal vector of 
the element. 
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At a point off the surface of the body, the components of the 
velocity are determined just as if the point of interest was a null point of 
a single element. The total velocity at the field point is the sum of the 
freestream velocity and the contributions of each of the elements of the 
body surface. The contribution of each of the elements is determined by 
calculating the influence coefficient based on the element geometry, and 
multiplying the result by the source density for the element. The total 
velocity at the field point may be expressed as 

N 

Vp = Voo + £ Cpq O q (125) 

q=i 



where p and q represent the field point and the source element 
respectively and the influence coefficient, 



C 



pq 





( 1 26) 



As discussed in section 6.0, the influence coefficient may be 
calculated by the exact method, or it may be approximated by the 
quadrupole or monopole method depending on the ratio of the distance, r 0 , 
between the field point and the centroid of the source element to the 
maximum dimension, t, of the source element. 
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The magnitude of the velocity at either the on-body or off-body 
points is given by 



|v 





(127) 



The pressure coefficient is calculated by using the result of 
equation (127) in equation (19), renumbered here for clarity. 



c p = 



P - P 



= 1 - 






CO 



(128) 



x» 
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11.0 STREAMLINE CALCULATIONS 



For the calculation of streamlines off the surface of the Dody, a 
timestep procedure is performed by calculating the velocity at the 
starting point of the streamline from equation (125), and advancing the 
streamline onetime increment by a fourth order Runge-Kutta integration 
to a new point (Ralston 1965). The timestep procedure is repeated, thus 
creating a streamline composedof finite segments. 

For the calculation of streamlines on the surface of the body, the 
streamline is started at a specified point and quadrilateral number. The 
local velocity is calculated from equation (124), and the values of a 
stream function are computed for each corner point. The stream function 
is chosensothat it has a value of zero at the last point on the streamline 
in the quadrilateral. The side of the quadrilateral through which the 
streamline exits is determined, and coordinates of the point on the side 
which has a stream function value of zero are computed. The direction of 
the streamline is verified by comparing it with the known direction of 
positive velocity. The next quadrilateral through which the streamline 
passes is determined by calculating the proximity of the new 
quadrilateral to the most recent point on the streamline. A circular area 
is computed which encloses the new quadrilateral with an additional 10% 
margin. If the last point of the streamline falls outside the circle, then 
- the quadrilateral is discarded and a new one selected until the streamline 
is adjacent to the new quadrilateral. 
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Figure 13 . Calculation of on-body streamlines 



This procedure is repeated along the surface of the body until all of 
the surface elements have been tested. The result is a streamline 
composed of segments from one side of an element to another. 
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12.0 DEVELOPMENT OF HIGHER ORDER PANEL METHODS 



The XYZ Potential Flow program assumes a constant element 
source panel as described in Section 3.3. Extensive use of the constant 
element source panel method has shown that the primary diasadvantage of 
the method is that, in order to obtain a highly accurate solution, a large 
number of surface elements must be used to discretize the body surface. 
The method has been applied to problems of increasingly complex 
configurations (Hess 1977). By doing so, the size of the coefficient 
matrix is increased resulting in increased computer time and cost. 
Additional cost is accrued due to the manhours required to prepare the 
input. Therefore, while the constant element methods have proven to be 
very successful, the cost has motivated the development of higher order 
methods. 

The higher order surface singularity methods discretize the body 
surface with curved elements having a variable source density, as 
compared to the flat elements of constant source strength used in the 
basic method. Hess (1973) showed that the effect of a curved surface and 
the effect of a variable source density are of the same order of 
magnitude. Therefore, the two effects must be used together to provide a 
"consistent” solution. The consistent higher order panel method provides 
the increased accuracy and speed desired for three dimensional Neumann 
- problems (Hess 1979). 

Accordingto Hess (1979), the evolution of the higher order panel 
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method from the constant element method involved the derivation of new 
influence coefficients based on the integration of a variable source 
density over a curved element. Other portions of the method were 
unchanged. However, the development of the higher order velocity 
equations also required different programming logic. 

In examining the potential for development of the higher order 
methods, Hess (1979) noted that "a consistent approach always uses a 
source polynomial one degree less than the panel polynomial." Through an 
independent effort, Brebbia (1984) presented a higher order approach 
using the direct method to solve for a surface potential polynomial 
stating that the potential function must be of a degree at least equal to 
the degree of the polynomial describing the element. Knowing that the 
velocity function is the derivative of the potential function, these two 
observations agree. Asa result of his derivations, Hess showed that the 
solution of a flat element with a constant source requires one integral, a 
paraboloidal panel with a linearly varying source density requires six 
integrals, and a cubic element with a quadratic source density requires 
twenty-three integrals. Development of higher order methods has 
focused on the paraboloidal element with the linearly varying surface, as 
solutions of higher order than that offered little benefit for the amount 
of effort required to produce a working program (Hess 1979). Hess 
(1979) and Eriksson (1983) have independently developed programs for 
three dimensional higher order panel methods. The higher order Hess 
program evolved from the constant element program which he developed 
in the early 1960s, while Eriksson developed a new program based on the 
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work of Johnson and Rubbert (1975). Continued work in the near future is 
expected to deal primarily with refinement of the paraboloidal element 
with a linearly varying source (Eriksson 1983). 

In order to alleviate the burden of preparing the input, a geometry 
package for input data generation has been developed which is 
incorporated into the Hess higher order panel program. This allows the 
user to enter relatively few points to describe the body. The geometry 
package enhances the surface representation by distributing additional 
points on the surface based on one of many algorithms or recurring 
geometries (Halsey 1978). 

As the state of the art in fluid dynamics has progressed, the XYZ 
Potential Flow program has seen increasingly complex applications 
requiring a great deal of effort in preparing the input, and requiring long 
computer run times. Hess (1979) reported the use of the Hess constant 
element program for a configuration utilizing 7000 effective elements. 
Realizing that the computation time increases as the square of the number 
of elements, it is easy to see the motivation for developing the higher 
order panel methods. Though modern computers offer storage capacities 
which can handle most applications of the constant element method, the 
higher order panel methods can provide equal accuracy for much less user 
effort. While the constant element method is still a versatile tool, 

- future generations of the surface singularity methods will be able to 
handle the more complex applications being demanded in fluid dynamics. 



71 



13.0 VFLOCITV CALCULATIONS FOR A TR1AX1AL ELLIPSOID 



As the only true Dody for which an analytical solution exists, a 
triaxial ellipsoid was selected for the sample calculations in order to 
compare calculated results with the analytical solution. Hess has made 
use of the triaxial ellipsoid throughout his works in developing both the 
constant element method and the higher order panel method. Therefore, 
the XYZPF program will be compared with existing results of the Hess 
method (Hess 1979). 

The triaxial ellipsoid utilized for the calculations has semiaxes 
dimensions of 1. 2, and 0.5 in the x. y, and z directions respectively. The 
surface was discretized by selecting fixed intervals of 0.1 in the y 
direction, and fourteen equal divisions of the 90° sector in the x-z plane. 
The values of x and z were then solved in terms of y and an angle 6. This 
method yielded 280 elements in the first octant fora total of 2240 
effective elements after employing symmetry. A FORTRAN program was 
used to generate the corner points and to prepare the input file for later 
use by the XYZPF program. 

Figures (16) and (17) show excellent correlation with the analytical 
solution and little difference from the Hess solution using 4320 effective 
elements. The use of the centroid as the control point is an 
- approximation used to simplify the multipole expansion of the potential 
about the origin of the local coordinate system. This approximation is 
valid for most elements. However, for elements which are long and thin, 
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the physical difference between the centroid location and the null point 
location is significant, and use of the centroid can produce significant 
error as may be observed in figure (16) when the value of y approaches 
2 . 0 . 



Recent calculations on the same body (Hess 1979) showed that 
results of at least equal accuracy could be obtained using only 480 
effective elements using the higher order panel method. These results are 
a significant demonstration of the value of the higher order panel method. 
Using the higher order panel method rather than the constant element 
method, the user has the option of obtaining equal accuracy with cruder 
discretization or higher accuracy for the same discretization effort. 

While the results of the triaxial ellipsoid show relatively little 
improvement in accuracy, the most significant advantages of the higher 
order panel method are evident fora body with concave regions (Hess 
1977). 
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ANALYTIC SOLUTION 

o HESS PROGRAM (4320 ELEMENTS) 
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Figure 16 . Comparison of analytic and calculated velocity distributions on an ellipsoid with 
axes ratios 1 :2:0.5. Velocities in the xz-plane. (from Hess and Smith 1962) 
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— ANALYTIC SOLUTION 
o HESS PROGRAM (4320 ELEMENTS) 




Figure 17 . Comparison of analytic and calculated velocity distributions on an ellipsoid with 
axes ratios 1 :2:0.5. Velocities in the yz-plane. (from Hess and Smith 1962) 



75 



14.0 CONCLUSION AND REMARKS 



The objectives of this paper were (1) to describe the details of the 
approximation of an arbitrary three-dimensional body using quadrilateral 
elements, and (2) to provide a detailed derivation of the exact source 
panel integrations. Both of these objectives were met. 

The method of surface discretization and source panel geometry is 
easily described using basic principles of geometry and vector algebra. 

By using quadrilateral surface elements, many surfaces can be discretized 
in a very straight forward logical fashion. The user can frequently 
visualize the contour lines of the surface which may be used to form the 
quadrilaterals, with some help from an intuitive approach to the fluid 
dynamics problem. The method of forming the planar quadrilateral 
element in the XYZPF Program differs slightly from the method presented 
by Hess and Smith (1962). The differences lie in the formation of the 
local coordinate system and the use of the centroid rather than the null 
point as the control point for applying the boundary conditions. The 
method of forming the local coordinate system has no effect on the 
potential flow calculations as long as one of the coordinate vectors is 
the outer normal to the planar element. The use of the centroid as the 
control point is an approximation used to simplify the multipole 
expansion of the potential about the origin of the local coordinate 
system. This approximation is valid for most elements. However, for 
elements which are long and thin, the physical difference between the 
centroid location and the null point location is significant, and use of the 
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centroid can produce significant error as may be observed in figure (16) 
when the value of y approaches 2.0. 

A detailed derivation of the exact source panel integration has not 
previously appeared in literature, though the results are summarized by 
Hess and Smith (1962). The derivations presented in this paper verify the 
equations presented by Hess and Smith (1962), and the equations used in 
the XYZPF Program. The integral expressions for the velocity 
components were evaluated exactly with no assumptions or 
approximations used in the course of the integrations. Since the method 
of integration reduces the surface integral to a line integral around each 
of the sides of the element, the integration method can be generalized for 
a planar element with any number of sides, though the surface 
discretization used in the XYZPF Program uses only quadrilateral 
elements. 

The calculation of potential flow about arbitrary three dimensional 
bodies is an engineering tool which is basic to design involving fluid 
dynamics. The XYZPF Program is a useful tool which has proven its value 
over the past 14 years. However, the increasing demands placed on this 
method are exposing the errors of the approximation as evident in the 
sample calculations presented in this paper. The requirement for 
increased accuracy has motivated the development of the higher order 
- panel methods. Some of the limitations imposed on the XYZPF Program 
were due to computer memory and speed limitations. Advances in 
computer performance may allow future investigators to eliminate some 
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of the simplifying approximations used in the XYZPF Program, allowing 
increased accuracy without violating computer limitations. Some 
modifications might include the use of the null point as the control point 
rather than the panel centroid (as is used in the Hess program), or 
extending the range in which the exact velocity calculations are 
performed. The gains in accuracy by modifying the "constant element 
method" are limited by the basic approximations of the planar element 
and the constant source density for each element. Significant gains are 
most evident in the higher-order panel methods. This author concurswith 
Eriksson (1983) in expecting advances in the surface singularity methods 
to focus on the "development and refinement" of the higher order panel 
methods. 
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oooooooooo 



APPENDIX I - XVZPF SECTION PF 1 



PROGRAM PFP 1 ( I NPUT= 1 2S , OUTPUT* 1 28 , TAPE5= I NPUT , TRPE6=0UTPUT , TRPE03 , 
1 TAPES-TAPE03 , TRPE04 , TRPE4=TRPE04 , TRPE50= 1 28 > 

XVZ POTENTIAL FLOW PROGRAM VERSION 4 SECTION 1 
READS INPUT AND COMPUTES QUADRILATERAL PARAMETERS 



FOR INFORMATION CONTACT 

BILL CHENG OR .JANET DEAN 
NUMERICAL FLUID DYNAMICS BRANCH CODE 1843 
NRUAL SHIP RESEARCH AND DEUELOPMENT CENTER 
BETHESDA, MARYLAND 20084 



C 



D I MENS I ON I NDEX(9, 3 > . 6(9,6 >, F(9 >, CZ(9 > , I P(9 >, XP(9 >, YP<9 ), 2P<9 ) 

1 , MSK( 100), WS<240>,PR0B(15),DM(650> 

COMMON X<800 >,Y(800 ),Z(800 ). 10(4 1,71 >, 8(250), T(4600),KP( 100) 
EQUIVALENCE <CZ< 1 >,F< 1 > > 

EQU I VALENCE < WS ( 1 ) , KP ( 1 ) ) , ( US ( 1 0 1 ) , MSK ( 1 )), (US(201 >, NP ) , ( US (202 ) , 

2 NSP > . < US < 203 ) , NEP > , ( US < 204 ) , NSE > . ( US ( 205 ) , M I X > , < US C 206 >,M! V > , 

3 < US ( 207 ) . M I Z ) , ( US (209 ) , I CUM). (US<210). ISM). (US(21 1 ).K), (US(213>, 

4 EPS), (US (208), IPS ). (US(2 12 ), IFF ), (US(217),XI ), (US(2 18 ) . VI ), 
5(WS(219),ZI ) 



EQU I VALENCE CV 12, V23 ) . ( V34 . Y4 1 ) 

I NTEGER P , P 1 , P2 . P3 . P4 PC , F'5 , P6 , F'7 . P8 , P9 
WRITE (6,5) 

5 FORMAT ( 4 YH 1 XVZ POTENTIAL FLOW PROGRAM SECTION 1, VERSION 4 ) 

10 FORMAT (II, 15R4 ) 

20 FORMAT ( IX. 1514) 

50 FORMAT ( IX, 2 1 7. 6E 12. 5) 

30 FORMAT ( IX, 5F 12.9) 



A. READ IN CONTROL PARAMETERS 
US (220 )=4 . 

K1=Q 



ID 1=0 



1 02=0 



I D3=0 
1 04=0 
1 05=0 



I D6=0 

107=0 

NAXN-70 

MAXM=40 

MAXNQE=650 

MAXPC=800 

i CNTRL= 1 

E0F50-0. 

READ (5 . 10)J,(PR0B(I ), 1=1, 15) 

IF (E0FC5) .EQ. 0. ) GO TO 9 
WRITEC6, 8) 

8 FORMAT (39H0N0 TITLE CARD FOUND - PROGRAM ABORTED ) 

STOP 

9 CONTINUE 
U=0 

WRITE (6,10) J,(PR0B(I ), 1=1, 15) 

SA= . 0 
SB=0. 

21 FORMAT (17H0N0. OF QUADS. =, 14 / 17H NO. OF SECTIONS® 

2I4/31H MAX. NO. OF ITERATIONS X FLOW , 13, 9H V FLOW , 
313, 9H Z FLOW, 13) 

READ (5 , 11)NQE,NSE,MIX,MIY,MIZ, ISM, EPS, IUCT, IPS. IFF, ISP 
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1 , I ED I T 1 , I ED I T3 , I ED I T4 , 1 TAPE , XCENTER , VCENTER , 2CEHTER 

11 F0RMRT(6I4.F8.5,8I4, 1X,3F5.3) 

IF (E0F(5) .EQ. 0. ) GO TO 10 
URITE(6, IS) 

18 F0RMRT(43H0N0 PARAMETER CARD FOUND - PROGRAM ABORTED ) 
STOP 



19 CONTINUE 

IF (IEDIT1.EQ. 1) I CNTRL=0 
l-JRITE (6,21 ) NQE . NSE, M I X, M I V , M 1 2 
31 FORMAT ("CONVERGENCE CRITERIR",F8.5) 

4 1 FORMAT ( 4H0 M , 7X , 2HX 1 , 1 2X , 2HX2 , 1 2X , 2HX3 , 12X , 2HX4 , 1 2X , 2HXP , 1 2X , 

1 2HXN, 12X, 1HR, 13X.3HCZ4/4H N,?X,2HV1, 12X,2HV2, 12X,2HV3, 12X,2HV4, 

2 12X , 2HVP, 12X , 2HVN . 12X , 2HFL , 12X . 3HC25/4H P, 7X , 2H2 1 , 12X, 2H22 , 12X , 

3 2HZ3 . 1 2X , 2HZ4 . 1 2X , 2HZF' . 1 2X , 2HZN , 1 2X , 4HCZ 1 , 1 0X , 3HCZ6 / ) 

42 FORMAT <4H 1 M, ?X, 2HX 1 , 12X, 2HX2 . 12X , 2HX3, 12X, 2HX4 , 12X , 2HXP , 12X, 

1 2HXN , 1 2X , 1 HR , 1 3X , 3HC24 /4H N, 7X, 2HY 1 , 1 2X . 2HY2 , 1 2X , 2HY3 , 1 2X , 2HY4 , 

2 12X, 2HVP , 12X . 2HVN, 12X, 2HFL , 12X, 3HC25/4H P, 7X, 2HZ 1 , 12X, 2HZ2, 12X . 

3 2H23 , 1 2X , 2H24 , 1 2X . 2H2P , 1 2X , 2H2H , 1 2X . 4HC2 1 ,1 OX , 3HCZ6 / > 

24 FORMAT <1H , I 1, 19H PLRNES OF SYMMETRY) 

240 WRITE (6,24) ISM 



270 WRITE <6,31) EPS 
280 IF (IPS.LE.O) GO TO 290 
285 WRITE (6,36) IPS, IF'F 

36 FORMAT (45H0NEM SOURCE DENSITY TO BE COMPUTED FOR QUADS., 1 4, 3H - , 
114) 

290 K=0 



WRITE(6,39)ISP 

39 FORMAT <9H0 ISP = , 13 ) 

HR I TE<6 . 37 > I ED I T 1 . I ED ITS, I ED IT4, I TAPE 
37 FORMAT (9H0IEDIT1 =. I3/9H I ED ITS =, I3/9H IEDIT4 =, I3/9H I TAPE =, 
1 13) 



WR I TE<6, 38 ) XCENTER , VCENTER , ZCENTER 
38 FORMAT < 1 OHOXCENTER = ' , F5 . 2/ 10H VCENTER = , F5 . 2 / 1 OH ZCENTER = , F5 . 2 ) 
MM=0 



MN=0 



P=1 



0 = 1.0 

DO 291 1=1.41 
DO 291 J-1.71 
291 IOC I . U)=0 
J=0 

C e. READ FIRST PT. 

I ERFi=0 

IF ( I TAPE . EQ . 1 ) GO TO 292 
2000 READ (5,40) X I , Yl ,2 1 , Nl .Ml , NS ,NE , UN 

IF <EOF(5).NE.O. .OR. NS.LE.O) GO TO 2050 
WRITE ( 6 , 45) I CNTRL , NS 
45 FORMAT (I 1 , 9H SECTION ,14) 

L I NE=0 
GO TO 293 

2050 IF< IERR.EO.O) GO TO 2200 
URITE<6,2100) 

2100 FORMAT (39H0N0 POINT CARDS FOUND - PROGRAM ABORTED ) 

STOP 

“ 2200 I ERR =1 
ITRPE= 1 

WRITE (6, 2300) 

2300 FORMAT (47H0ERR0R IN INPUT - POINT CARDS NOT ON INPUT FILE , 10X, 

1 53HPRQGRRM WILL CHANGE I TAPE TO 1 AND TRY TO READ TAPE50 ) 

IF ( EOF (5 ) . EQ . 0 ) UF;ITE(6,2400) XI, VI. 21 

2400 FORMAT ( 1 1H0EXTRA FLOW, 3F 12.5. 5X.20HH ILL NOT BE COMPUTED ) 
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292 RERD(5G, 40) XI , V I . 2 1 , Nl , hi , NS . NE , UN 
40 FORMAT (3F 12. 9, 4 1 4, F 12. 9) 

E0F50=E0F(50) 

IF <E0F50.ME.0. .OR. MS.LE.O) GO TO 2450 
WRITE <6, 45) ICNTRL.NS 
L I NE=0 
GO TO 293 

2450 IF < I ERR. EC!. 0) GO TO 2500 
WRITE(6,21Q0) 

STOP 

2500 I ERR= 1 
I TRPE=0 
WRITE <6, 2600) 

2600 FORMAT ("ERROR IN INPUT - POINT CAROS NOT ON TAPE50", 10X, 

1 "PROGRAM WILL CHANGE I TAPE TO 0 AND TRV TO READ INPUT FILE") 

GO TO 2000 

293 UNR-UN 
NSS-NS 
F'C= 1 

IF <NE . EQ . 0 ) GO TO 2700 
MM I N=N I 
MMAX=N I 
NM I N=M I 
NMRX=M I 
GO TO 300 
2700 MM IN-MI 
MMAX=M I 
NM I N-N I 
NMRX=N I 
GO TO 300 

295 I F < I TAPE . EO . 1 ) GO TO 297 

READ (5,40) XI, VI, 21, Nl, Ml, NS, ME, UN 
IF (E0F(5 ).EQ .0. ) GO TO 299 
NS=0 
X I =0 . 

V I =0 . 

21 = 0 . 

GO TO 299 

297 READ (50, 40 ) X I , V 1 , 2 1 , N I . M I . NS . ME, UN 
E0F50=E0F(50) 

IF (E0F50 . NE . 0 ) NS=0 

299 PC=PC+ 1 

IF (NS.NE.NSS) GO TO 330 

300 IF (NE.EQ.O) GO TO 304 

301 I W=N I 
N I =M I 
M I = I W 

C C. STORE PT. IN PT. ARRAY 

304 IF (MRXPC+ 1— PC ) 295,305,310 

305 WRITE (6, 306) NS, Ml, N I 

306 FORMAT (60H ERROR IN INPUT - THERE ARE TOO MANY DATA POINTS IN 
IT ION , I4,30H - POINTS BEGINNING WITH M =, 14, 5H N = , 14, 

2 17H WILL BE IGNORED ) 

LINE-LINE+1 
I D4= I D4+ 1 
GO TO 295 
310 X(PC)=XI 
Y(PC >=Y I 
2(PC)=2I 

IF (MI.LE.MAXM .AND. NI.LE.MAXN) GO TO 315 
WF;ITE(6.31 1 ) MI,NI 



85 



L I NE=L I NE+ 1 

311 FORMAT OSH ERROR IN INPUT - I HURL 1 0 M, N INDICES , 10X, 

1 14HP0INT WITH M = , I4,5H N = , 14, 17H WILL BE IGNORED ) 

I D5= I D5+ 1 
PC=FC- 1 
GO TO 205 
315 ID(11I,NI>=PC 

MMRX=11RX0(MMRX.MI ) 

MMIN=11IN0<1111IN,11I ) 

NMRX-MRXO<NMRX,HI > 

Nil 1 11=11 1 NO (Nil IN, HI ) 

GO TO 295 

330 IF (IEDIT1.EO. 1) GO TO 294 
IF (LINE.LT.40) GO TO 333 
WRITE(6,42> 

L I NE=0 
GO TO 294 

333 URITE<6,41 ) 

294 CONTINUE 

C E. DO LOOPS TO SWEEP PT. RFiRflV 

N 1 =Nf1 1 N 
MM2=Mf1flX-MM I N 
NN2=NMRX-NM I N 

IF < MOD(11M2,2).EQ.O .AND. M0D(NN2.2>.EQ.0> GO TO 332 

WR I TE (6 , 33 1 > NSS . MM I N . 11I1AX , NM I N . NMRX 
L I NE=L I 11E+ 1 

331 FORMAT ( 16H0ERR0R - SECTION , I5,45H DOES NOT HRUE QUADS ARRANGED IN 
1 BLOCKS OF 4 , 9H MHIN= ,12. 6H HMAX= , I2,6H NM I N= . 12 , 

25H NMAX= ,12) 

106=106+1 

332 MM2-MH2/2 
11112=11112 /2 

DO 404 NN= 1 , HN2 
HI -MM IN 

DO 402 MM= 1 , MM2 
110= 1 

C F. HRUE 9 CORNER PTS. BEEN G I OEM 

I T= I D(M 1 , N 1 >* I D(M 1+ 1 . N 1 >+ 1 0<M 1+2 . N 1 >+ 1 D(M 1 . N 1+ 1 )* 1 0(11 1+ 1 , N 1+ 1 >+ 

1 I D <11 1 + 1 , N 1 +2 )+ 1 D < M 1 . N 1 +2 >* 1 0 ( M 1 + 1 . N 1 +2 >+ 1 D <11 1 +2 , N 1 +2 > 

IF( IT.EQ.O ) GO TO 402 
I EFiR=0 
112=111+1 

N2=N 1 + 1 =111,112 
DO 400 N=N 1 , N2 
GO TO <334,335.336,33?) NO 

334 P 1=1 0(11 , N ) 

P2=ID<M+1,N ) 

P3= ID (11+1,11+1) 

P4=ID(M . N+ 1 ) 

P5= 1 0(11+2,11 ) 

F'6= 1 0(11+2, N+1) 

P7= ID (11+1,11+2) 

F'S=ID(H , N+2 ) 

P9=F' 1 

I F < <X<P 1 ) . NE . X(P2 ) . OR . V<P 1 ) . HE V(P2 > . OR . Z(P 1 ) . NE . 2<P2 ) ) . AND . 

1 (X(P1 ).NE.X(P4).0R.V(P1 ).NE.V(P4).0R.Z(P1 ).NE.Z(P4))> GO TO 340 

F'9= 1 0(11+2, N+2) 

GO TO 340 

335 P 1 = 1 0(11 ,11+1) 

P2= ID(11 ,N ) 

P3=ID(M+1 ,N ) 
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P4=IDCM+1 J N+1> 

P5= IDCM , H- 1 > 

P6=IDCM+1,N-1) 

P7=IDCM+2,N ) 

PS=!DCM+2.N+1 ) 

P9=P 1 

I F< <X<P 1 > . NE . XCP2 ) . OR . V<P 1 ) . ME VCP2 > OR . ZCP 1 > . ME . ZCP2 ) ) RMD . 

1 CXCP 1 ). NE.XCP4) .OR. VCP1 ) . NE . VCP4 ) . OF; . 2 < P 1 ). NE . ZCF'4 )>) GO TO 340 

PQ=I0CH+2,N-1> 

GO TO 340 

336 P1=IOCf1+1,N ) 

P2= I DCM+ 1,N+ 1 ) 

P3=IDCM ,N+ 1 ) 

P4=ID(H jN ) 

P5=IO<M+1,N+2> 

P6=ID<M ,N+2) 

P7=IDCf1-1,N+1) 

P8=IQCM-1 ,N > 

P9=P 1 

I FCCXCP 1 ) . ME . X<P2 ) . OR . VCP 1 > . NE . VCP2 ) . OR . ZCP 1 ) . HE . ZCP2 ) > . RMD . 

1 C X C P 1 ) . ME . X C F'4 ) . OR . S' C PI). HE . V C P4 ) . OF; . Z < P 1) . NE . Z C F'4 ) ) ) GO TO 340 

P9=l0CH-1,H+2> 

GO TO 340 

337 P1=IDCM+1,N+1) 

P2-IDCM M+ 1 ) 

P3= IDCM ,H > 

P4=IDCH+1,N ) 

P5= 1 0<t1- 1 ,H+ 1 ) 

P6=I0CM-1 ,N > 

P7=ID(f1 j H- 1 ) 

P8= I D<M+ 1 , N- 1 > 

P9=F 1 

I F C C X C P 1 ) . HE . X C P2 ) . OR . V C PI). HE . V C P2 ) . OR . 2 C P 1) . HE . Z < P2 ) ) . HMD . 

1 CXCP 1 > . HE . X<P4 > . OR . YCP 1 > . HE . YCP4 ) . OR . ZCP 1 > . HE . Z(P4 > > > GO TO 340 

P9= 1 0 < M— 1 , H- 1 > 

340 IP<1)=P1 
I P ( 2 )=P2 
I PCS )=P3 
IPC4)=P4 
I PC 5 )=P5 
IPC6)=P6 
IPC7>=P7 
I PCS >=P8 
IPC9)=P9 

G2 COMPUTE NORMAL UECTOR CXH.VH.ZM) 

X1=XCP3)-XCP1 ) 

X2=XCP4)-XCP2) 

V1=VCP3)-VCP1 ) 

V2=VCP4 )-VCP2 ) 

Z1=ZCP3)-ZCP1 ) 

Z2=ZCP4 )-ZCP2 ) 

XN=V2*2 1-V 1+Z2 
VH=X 1+Z2-.X2+Z 1 
ZN=X2*V1-X1*Y2 
R=SQ2CXH,VN,ZH) 

IF CR.GT. .00000000001) GO TO 345 
UR I TEC 6, 343) 

343 FORMAT C33H ERROR IH INPUT - ZERO AREA QUAD ) 

L I HE=L I NE+ 1 
ID1= ID1+ 1 
AQ=0 . 
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XC-O. 

YC=0. 

ZC:=0 . 

FL=0. 

cza>=o. 

CZ(4 >=0 . 

CZ<5)=0. 

CZ(6)=0. 

I ERR= 1 
GO TO 351 
345 CONTINUE 
XN=XN/R 
VN=YN/R 
ZN=ZN/R 
RO= . 5*R 

C COMPUTE CENTRO I D 

X1=X<P3 )-X(P2) 

V1=Y(P3>-Y<P2> 

Z1=Z<P3)-Z<P2> 

X5=V 1 +Z2-Y2+Z 1 
V5=Z 1+X2-Z2+X 1 
Z5=X 1*V2-X2 :+: Y 1 
ft 1 =SC|2 < X5 , Y5 , Z5 ) 

R2=F;-fi 1 
I T= 1 

XC=(X(P2 )+X<P4 >+<R1*X<P3 )+R2*X( P 1 > )/R >/3 . 

VC= < V < P2 )+ V < P4 >+ < ft 1 +V < P$ >+R2+Y ( ? 1 > ) /R > /3 . 
ZC=(Z<P2 )+Z(P4 >+<R 1+ZCP3 HA2+Z<P 1 ) )/R )/3. 

C COMPUTE SECOND RND THIRD UECTORS 

045 X4=VN+Z1-V1+ZN 
Y4=ZN+X 1-Z 1*XN 
Z4=XN*V 1-X 1*VN 
fl= 1 . /SQ2CX4 . Y4 , Z4 ) 

X4=X4*A 

Y4=V4*H 

Z4=Z4+R 

X3=ZN*V4-Z4*VN 

Y3=XN+Z4-X4+ZN 

Z3=VN*X4-V4*XN 

C COMPUTE POINTS IN QUAD SYSTEM 

DO 947 1=1,9 
L= I P< I > 

XP< I >=X3*(X(L )-XC: >+V3*(V<L )-VC >+Z3+(Z(L >-ZC > 
VPC I )=X4'*'CX(L >-XC >+V4*(V<L )-VC )+Z4 + (Z(L )-ZC ) 
94? ZP < I >=XH+ < X ( L >-XC )+YN* ( V < L )- VC >+Zti+ (ZCL )-ZC > 
C COMPUTE MATRIX COEF. TO FIND SURFACE EQ. 

DO 949 1=2,9 
G( 1 , 1 )= 1 . 

G( I ,2)=XP< I ) 

G( 1 ,3>=VP( I ) 

0(1,4 )=XP( 1 >**2 
G( 1 ,5>=VF’( I )**2 
G< I ,6)=VP(I >+XP ( I ) 

949 F( I )=ZP( I ) 

DO 953 1=1,6 
G( 1 , I >=GC9, I ) 

G<5, I >=G<5, I )+G(6, I > 

953 G<6, I >=G<7, I >+G<8, I > 

F( 1 )=F(9> 

F(5)=F(5)+F(6) 

F(6)=F(?)+F<8) 
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C SClLUE MATRIX EQ. G+CZ=F FOR CZ 

CALL MATINS(G,9,6,F,6, 1, DETERM, I DM, INDEX) 

IF (IDN.EQ. 1) GO TO (955,900) IT 
I ERR= 1 

WRITE (6, 954) 

954 FORMAT (33H ERROR IN INPUT - SINGULAR MATRIX 
LINE=LINE+1 

102 = 102+1 
GO TO 960 

955 IT=2 

C FIND NEW NORMAL UECTOR 

XN=XN-CZ<2 )+X3-CZ<3 )*X4 
VN=VN-CZ(2 )*V3-CZ<3 )*V4 
ZN=ZN-C:ZC2 )+Z3-CZ<3 )*24 
A=1./SQ2(XN,VN,ZN) 

XN=XN+fl 
YN=VN+A 
2N=ZN'+A 
GO TO 945 

C STORE DATA 

960 B(J+1)=XP(1) 

S<J+2)=VP(1) 

B(J+3)=XP(2) 

BCJ+4 )=YF'(2 ) 

B(J+5 )=XP(3) 

BCJ+6 )=XP(4) 

B(J+? )=VP<4 ) 

B( J+S >=X3 
B( J+9 )=V3 
B(J+ 10 )=Z3 
B(J+ 1 1 >=X4 
B( J+ 12 )=Y4 
B(J+ 12 )=Z4 
B< J+ 14 )=CZ( 1 ) 

B(J+15>=C2<4> 

B(J+ 16)=C2<5) 

B(J+17)=CZ(6) 

IF (K .LT. 7+MRXNQE ) GO TO 965 

107=107+1 
K 1=K+K 1 
K=0 

965 CONTINUE 
T(K+ 1 )=XC 
T(K+2 )=VC 
T (K+3 )=2C: 

T(K+4 >=XN 
T(K+5 )=YN 
T(K+6)=ZN 
T(K+7)=AQ 

C COMPUTE QUADRUPOLE MOMENTS 

XI 1=XP< 1 )+XP<2) 

XI2=XF'< 1 )+XP(4 ) 

X 1 3=XP < 3 >+XP < 2 ) 

XI4=XP(3)+XP(4 ) 

XI5=XP(2)+XF’<4) 

VI 1=YP( 1 )+VP(2 ) 

YI2=YP< 1 )+YP(4 ) 

YI3=YP(3)+YP<2) 

YI4=VP(3)+YP(4) 

V 1 5=VP<2 )+YP(4 ) 

Fl1=A1/24. 
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R2=fl2/24 . 

R3=RQ/12. 

RXX= (XII **2+X 1 2**2 >*R 1 + ( X 1 3**2+ X 1 4**2 >*R2+X 1 5**2*R3 
RXV= ( X II *Y II +X 1 2*Y 1 2 )*R 1 + ( X I 3*V 1 3+X 1 4*V 1 4 >*R2+X 1 5*V 1 5*R3 
RYY=< Y I 1**2+ Y 1 2**2 >*R 1+(Y 1 3**2+ Y 1 4**2 )*R2+V 1 5**2*R3 
C COMPUTE SOLID ANGLE 

XX=XC-XCENTER 
YV=YC-VCEMTER 
ZZ=ZC-ZCENTER 
X 1 =XX*X3+ YY+Y3+Z2+23 
V 1 =XX*X'4+ VV* V4+ZZ*Z4 
2 1=XX*XN+VY*YN+ZZ*ZN 
F;D= 1 . /302(X 1 , V 1 , 2 1 ) 

RCU=RD**3 

R8U=RCU++'2*RD 

Sfl=SR+Z 1*(RO*RC.U-( CfiXX+XY 1**2+Z 1++2-4 . *X 1**2 > 

1 +RYY* < X 1 **2+2 1 **2-4 . *V 1 **2 ) )* 1 . 5- 1 5 . *X 1 *V 1 *fiXV >*RSU ) 

B(.J+ 18 >=RXX 
B(J+19>=flXY 
B( J+20 )=RVV 
C ERROR TESTS 

D 1 =SQ2 < ( XP ( 3 >-XP < 1 ) ) , < VP < 3 )-VP ( 1 ) ) , 0 . ) 

D2=SQ2 < < XP ( 4 :>-XP (2)), ( YP < 4 >-VP < 2 > ) , 0 . ) 

FL=.5*RNRX1(01,D2) 

0223=063 ( CZ(2) ) + RB8C C2(3> ) 

IF( RBS(CZ<2 ) )+RB8(CZ<3 > ) .GT. FL*.001 ) 00 TO 970 

IF < RBSC CZ( 1 )) .LT. FL+.3 ) GO TO 977 
970 UR I TE <6 ,975) CZ23 

975 FOFlMRT (29H QUESTIONABLE POINT -POOR FIT .6E14.3) 

I EFlFl= 1 
L IHE=LINE+ 1 

977 IF (VP( 1 )*XP(2 )-YP(2 >*XP( 1 ) .GE. 0. .AND. 

1 YP C 2 )*XP ( 3 )-YP C 3 >*XP ( 2 ) .GE. 0. .AND. 

2 YP<3 )*XP<4 )-YP(4 >*XP<3 > .GE. G. .AND. 

3 VP ( 4 )*XP < 1 >-VP ( 1 )*XF' ( 4 > . GE . 0 . ) GO TO 984 

980 UF!ITE(6, 1000) ( XP< I ), VP< I ), 1 = 1,4) 

1000 FORMAT (4 1H ERROR IN INPUT - CROSSED OR CONCRUE OUftD 
1 4(2F 10.5.3X) ) 

I ERR= 1 
L I NE=L I NE+ 1 
I D3= I D3+ 1 



984 CRCF=SQ2 < < XP < 2 >-XP < 1 ) ) , ( YP < 2 >-YP < 1 > ) , 0 )+XP < 3 )-XP < 2 )+ S02 < < XP < 1 >- 
1 XP(4>), <VP< 1 )-YF'<4)),0. )+SQ2( (XP(4 )-XP(3 ) ), (VP(4)-VP(3)),0. ) 

IF < 36 . *RQ .GT. CRCF**2 ) GO TO 986 
L I NE=L I NE+ 1 
UR I TEC 6, 96 1 ) 

961 FORMAT (24H UfiRNING LONG THIN QUAD.) 

986 IF ( 21 .GE. 0. ) GO TO 351 

350 WRITE (6,85) 

85 FORMAT (35H QUESTIONABLE POINT - INWARD NORMAL) 

I EFlR= 1 
L I NE=L I NE+ 1 

C J. EDITE OUftD INFORMATION 

351 IF ( I ED I T 1 . EQ . 2 . AND . IERFi.EQ.0) GO TO 354 
IF ( IEDIT1 . EQ. 1 ) GO TO 354 

GO TO (356.357,358,359) NQ 

356 MF;ITE(6,51 ) M,X(P1).X(P2),X(P3),X(P4).XC.XN,RQ ,CZ(4) . 

1 H . Y(P 1 >, Y(P2 ), Y(P3 > , Y(P4 ), YC, YN , FL ,CZ(5) ’ 

2 P, Z(P 1 ), Z(P2 ), Z(P3 )', Z(F'4 ), ZC, ZN‘ CZ( 1 ), CZ(6 ) 

GO TO 360 

357 WRITE (6,51) M,X(P2),X(P3),X(P4),X(P1 ).XC.XN,flQ ,CZ(4) , 
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2 

GO TO 360 




358 WRITE<6,51 > 




2 

GO TO 360 



359 URITEC6, 5 1 ) 




360 CONTINUE 
51 FORMAT <1H 




L I NE=L I NE+4 

IF <LINE.LT. 50) GO TO 354 
352 WRITE (6,42) 

L I NE=0 

354 CONTINUE 
J=J+20 
!=P 

DM< I >=UNR 
F-F+ 1 
NQ=NQ+ 1 
I ERR=0 
349 K=K+7 

C K. WRITE OUT BLOCK OF B RRRRV IF FULL 

IF (U.LT.240) GO TO 400 

355 WRITE (04) Q, <B< I >, 1 = 1,240) 



400 CONTINUE 
402 M 1=M 1+2 
404 N1=N1+2 

11. SET FOR NEXT SECTION 



DO 405 N=NM ! N, NfIRX 
405 IDai.N)=0 
PC= 1 

IF (ME . EO . 0 ) GO TO 410 
MMRX=NI 

nmN=m 

NM I N=N I 
NMRX=M I 
GO TO 420 
410 MMAX=MI 
MM I N=N I 
NM I N=N I 
NMflX=N I 
420 NE=ME 
UNR=UN 

IF (NS.LE.O) GO TO 500 
WRITE<6,45> I CNTFiL , NS 
L I NE=0 
GO TO 300 

500 WRITE (04 > Q, <B< I ), 1 = 1,240) 
550 NP=(K+K1 )/7 
I SMP = ISM + 1 

GO TO <595,590, 580 , 570 > , I SMP 
570 Sft=Sft+SA 
580 SR=Sft+SR 




C L. END OF DO LOOP GUER PT. RRRRV 
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590 SA=SR+sA 

C 01 WRITE PARAMETERS AMD T ARRflV OH TAPE 

595 J = 1 

IF < ITAPE.EQ. 1 .AMD. E0F50 .HE.O ) GO TO 601 

597 WS(J) = XI 
WS(J+20) = Vi 
WS(J+40> = Zl 
J = J+1 

IF(XI**2 + Vl**2 + Zl**2) 599,599.598 

598 URITE(6,600) XI ,VI ,ZI 

600 FORMAT ( 1 1H0EXTRA FLOW, 10X,3F 12.5) 

601 READ(5,4Q) XI .VI .21 

IF (E0F(5> EQ. 0) GO TO 597 

X 1 =c« . 

V I =0 . 



21 = 0 . 

GO TO 597 

599 IF (ISP.LT.O) GO TO 605 

WRITE (03) <PROB<! ), 1=1, 15) 

WRITE <03 ) <WS(I ), 1 = 1,220). I ED ITS. IEDIT4 
WRITE (03) (T(l), 1 = 1,10 
WRITE (03) (011(1 ), 1 = 1. HP) 

605 CONTINUE 

C HI CHECK SOLID ANGLE 

610 WRITE (6.80) 3A 
80 FORMAT ( 14H0S0L ID ANGLE = ,F8.3) 

620 REWIND 04 
REWIND 03 



622 I D3= I D 1 + 1 D2+ 1 03+ 1 04+ 1 D5+ 1 D6+ 1 D7 

IF (IDS.EQ.O .AND. NP.EQ.NQE) GO TO 638 
WRITE(6,625) 

IF ( ID1 . GT.O) WRITE (6, 628) ID1 
IF (ID2 .GT.O) WRITE(6,629) ID2 
IF (IDS .GT.O) WRITE (6. 630) IDS 
IF ( ID4 .GT. 0) WR I TE(6, 63 1 ) ID4 

IF (IDS .GT. 0) WRITE (6, 632) IDS 

IF (!D6 . GT. 0) WRITE (6, 633) 106 

IF (1 07 .GT. 0) WRITE (6, 634) NP, MflXNQE 

IF (NP.NE.NQE) WRITE (6. 637) HP,NQE 
STOP 



625 FORMAT (38H0FATAL 

628 FORM AT (1 HO, I5,31H 

629 FORMRTOHO, I5,44H 

630 FORHAT< 1 HO, I5/26H 

631 FORMAT ( 1 HO, 15, 32H 

632 FORMAT (1H0‘ I5,34H 

633 FORMRTOHO. I5,52H 

IF 4 ) 

634 FORMRTOHO, I5.48H 

1 14) 

637 FORMAT OHO, I5.28H 

638 IF (ISF.LE.O) GO 



ERROR IN DATA - PROGRAM ABORTED) 

QUADRILATERALS WITH ZERO AREA > 

QUADRILATERALS GENERATE A SINGULAR MATRIX 
CROSSED QURDFl I LATERALS ) 

SECTIONS HAVE TOO MANS-' POINTS ) 

POINTS HRUE I HURL ID M.N INDICES ) 

SECTIONS DO NOT HRUE QUADS ARRANGED IN GROUPS 

QUADRILATERALS GiUEN, EXCEEDING THE LIMIT OF 

QUfiDR I LATERALS G I UEN , NOT ,14) 

TO 640 



WRITE (6. 639) ISP 

639 FORMAT (7H0 ISF- , 14, 19H - PROGRAM ABORTED ) 
STOP 

640 CONTINUE 

C 02 READ PEPS2 AND TRANSFER TO IT 

STOP 1 
END 



C 



FUNCTION SQ2(X,V,Z> 



> 



0 
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C COMPUTE SQUflR ROOT OF R**2 

R= ABSCX >+ABS(V ) +ABSCZ) +.0000000000001 
700 RS=X**2+V**2 +2**2 
R=R+RS/R 
F;= . 25+R+RS/R 
R=R+RS'/R 
S02=.25*R+RS/R 
RETURN 
END 

SUBROUTINE MfiT I NS (A , NR , N 1 , B , NC , M 1 , DETERM , ID, INDEX) 

MATRIX INVERSION UITH RCCOMPRNVING SOLUTION OF SIMUL. EQ. 

PI MOT METHOD 

FORTRAN IU SINGLE PRECISION WITH ADJUSTABLE DIMENSION 
NOVEMBER 1971 S GOOD NRURL SHIP R & D CENTER 
WHERE CALLING PROGRAM MUST INCLUDE 

DIMENSION A(NF; , NR > , B(NR.NC), I NDEX(NFl , 3 ) 

WHERE NR , NC ARE DIMENSIONS OP A, B, INDEX 
N 1 IS THE ORDER OF A 

Ml IS THE NUMBER OF COLUMN VECTORS IN B ( MR',' BE 0> 
DETERM WILL CONTAIN DETERMINANT ON EXIT 
ID WILL BE SET BV ROUTINE TO 2 IF MATRIX A IS 
SINGULAR, 1 IF INVERSION WAS SUCCESSFUL 
MRTRIX A' (INPUT MATRIX) WILL BE REPLACED BV A INU 
MATRIX B: THE COLUMN VECTORS WILL BE REPLACED 
BV CORRESPONDING SOLUTION VECTORS 
INDEX: WORKING STORAGE ARRAY' 

IF IT IS DESIRED TO SCALE, THE DETERMINANT CARD 29 MAV BE 
DELETED AND DETERM PRESET' BEFORE ENTERING THE ROUTINE 

DIMENSION A (NR, NR), B(NR.NC), INDEXCNR.3) 

EQU I VALENCE ( I ROW , JROW ) , ( I COLUM , JCOLUM ) , ’ ( RMflX , T , SURF ) 

INITIALIZATION 

N=H 1 

n=n i 

DETERM= 1 . 0 
DO 20 J= 1 , N 

20 INDEXCJ,3)=0 
DO 550 I = 1 , H 

SEARCH FOR PIVOT ELEMENT 

AMfiX =0.0 
DO 105 J= 1 , M 

I F < !N0EX(J,3)-1 ) 60, 105, 60 
60 DO 100 K= 1 , H 

IF( INDEX(K,3)- 1 ) 80, 100, 715 
SO IF ( AMfiX -ABS (A(J,K))) 85, 100, 100 
85 I ROW = J 
I COLUM = K 
AMfiX = ABS (R(J,K)) 

100 CONTINUE 
105 CONTINUE 

INDEX Cl COLUM. 3) = I NDEXt I COLUM . 3 ) + 1 
INDEX (I, 1) = I ROW 
INDEX (I, 2) = I COLUM 

INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 



o o o o no o o o 



IF C I ROW- 1 CQLIIM > 140, 310, 140 
140 0ETERM= -DETERM 
DO 200 L=1,N 
SURF- ACIROU.L) 

R < I ROW , L )=R C I COLUM , L > 

200 RC I COLUM , L )=SURF' 

I F < M > 310, 310, 210 
210 DO 250 L= 1 , M 
SURP=BC IROU.L) 

B C I ROW , L >=B C I COLUM , L ) 

250 B< I COLUM, L )=SURP 

DIO IDE PIUOT ROW BV PIUOT ELEMENT 

310 PIOOT = RC I COLUM, I COLUM) 

DETEFlM=DETERM+P I UOT 
330 fl< I COLUM, I COLUM) = 1.0 
DO 350 L= 1 , M 

350 R ■: I COLUM ,L)=RC I COLUM , L ) /P I UOT 
IF CM) 380, 380, 360 
360 DO 370 L=1,M 

370 B< I COLUM , L )=B < I COLUM , L > /P I UOT 

REDUCE NOH-PIUOT ROUS 

3S0 DO 550 L1=1,N 

I FCL 1-1 COLUM) 400, 550, 400 
400 T=RCL1 , I COLUM) 

RCL 1 , i COLUM )=0 . 0 
DO 450 L= 1 , N 

450 RCL 1 , L )=RCL 1 , L )-RC I COLUM, L )*T 
IF CM) 550, 550, 460 
460 DO 500 L= 1 , M 

500 BCL 1 , L )=BCL 1 , L >-B< I COLUM, L >*T 
550 CONT I HUE 

INTERCHANGE COLUMNS 

DO 710 1=1, N 
L=N+ 1-1 

I F C I NDEX ( L , 1 >- INDEX < L , 2 > ) 630 , 710, 630 
630 JROM= INDEXCL, 1) 

JC0LUM= I NDEX CL . 2 ) 

DO 705 K= 1 , N 
SURF' = RCK'JROU) 

R C K , JROU )=fl C K , JCOLUM > 

RCK , JCOLUM )=SURP 
705 CONTINUE 
710 CONTINUE 
DO 730 K= 1 , N 

IFC INDEXCK,3)-1 ) 715, 720, 715 
720 CONTINUE 
730 CONTINUE 
ID= 1 

810 RETURN 
715 1 0=2 

GO TO 810 
END 
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APPEND I X II - IVZPF SECT I OH PF2 



PROGRAM PF2(0UTPUT= 128, TRPE6=0UTPUT, TRPE03, TRPE02, TRPE08, 

1 TAPEOQ, TflPE0 1 , TAPE 1 1 , TAPE04, TRPE4=TAPE04 , 

3 TRPE1=TRPE01, 

2 TRPE3=TRPE03 , TAPE2=TAPE02, TAF'E8=TAPE0S, TRPE9=TRPE09 ) 

XV2 POTENTIAL FLOW PROGRAM VERSION 4 SECTION 2 
COMPUTES MATRIX COEFFICIENTS 



C 



C 



COMMON EX 24 1 ), T(460Q >, U 1 ( 1000 ), V2C 1000 >, V3( 1000 ), C 1 ( 900), C2( 9 
200), CSC 900), US C300 ) , PFiOB CIS), 

3VX < 8 > , V V ( 8 ) , 02 C 8 ) 

EOU I VALENCECY3, Y2 > , (US (20 1 ) . MP ) . (USC2 10 ), SVM ), (US(2 1 1 >, KM ) 

1 , ( US < 2 1 2 > , I PF ) , ( US ( 208 > . I PS > , ( KM . NK > 



INTEGER SVM 
BLK= 1 . 0 



IDU=0 

READ (03) (PROBd ), 1 = 1, 15) 

URI TEC6,5> 

5 FORMAT (49H0XVZ POTENTIAL FLOW PROGRAM SECTION 2. VERSION 4 ) 

90 FORMAT ( 1H0, 15A4) 

UR I TE (6 , 90 ) CPR'OBC 0,1 = 1,15) 

A. READ PARAMETERS, T ARRAY, FIRST BLOCK OF B ARRAY 

READ (OS') (US (I ), 1 = 1,220) 

READ (03) <T< I >, 1=1 ,MK> 

READ (04) (B(i ), 1 = 1,241) 

B. START LOOP OVER QUAORALATERRLS 
K= 1 



•J=1 

P=1 



JC= 1 



JV=2 



NPN=5*< CNP+4 )/5 ) 

KMM=7+MPN+ 1 

290 IFCBC 1 )-P >595, 295, 595 

98 FORMAT ( " PO I NTS OUT OF ORDER B( 1 )=“ , 1F4 . 0 , " F'=" , 1F4 . 0 ) 

295 J=2 

296 X 1=B(J) 



V1=B(J+ 1 ) 
X2=B( J+2 ) 
Y2=B< J+3) 
X3=B(J+4) 



X4=B( J+5) 

V4=B(J+6) 

XN=T(K+3) 

YN=T (K+4 ) 

2N=T (K+5 ) 

XC=T(K> 

VC:=T(K+ 1 ) 

ZC=T(K+2) 

A=T (K+6) 

XX=B( J+7 ) 

VX-BCJ+S) 

ZX=B(J+9 ) 

XV=B(.J+ 10 ) 

YY=B(.J+11) 

ZV=B(J+12) 

C Cl COMPUTE LENGTH OF SIDES OF QUAD 

0 12=SQ2F(X 1 , X2, V 1 , Y2 . 0 . , 0 . ) 
D23=SQ2F(X2 . X3 . Y2 . Y3, . 0 . . 0 ) 
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034=S02F<X3,X4,Y3,Y4. .0, .0) 

041=SQ2F<X4,X1,Y4,Y1, .0, .0) 

C C2 COMPUTE SLOPE OF SIDES 

I F < X2-X3 >305 , 300 , 305 
300 0123=1. 

GO TO 310 

305 CM23=<Y2-Y3 >/(X2-X3 ) 

0123=0. 

310 IF<X3-X4>315, 31 1,315 

311 0134=1. 

GO TO 320 

315 CM34= < V4-Y3 )/<X4-X3 ) 

0134=0. 

320 I F(X4-X 1)325,321,325 

321 0141=1. 

GO TO 330 

325 CM41=CV1-V4)/(X1-X4> 

0141 = 0 . 

330 I F < X 1 -X2 >335 ,331, 335 

331 Cl 12=1. 

GO TO 340 

335 CM 12=CV2-V 1 >/<X2-X 1 > 

Cl 12=0. 

C 03 COMPUTE QURDRRPOLE MOMENTS 

340 Cl XX=B ( J+ 1 7 
CIXY=B<J+18> 

CIVV=B<J+19> 

CV 12=0.0 
0X12=0.0 
CY23=0 . 0 
0X23=0 . 0 
CV34=0 . 0 
CX34=0 . 0 
0Y4 1=0 . 0 
0X41=0.0 

0 04 COMPUTE SIM ft NO COS OF SLOPE ANGLE FOF; EACH SIDE 

IF < D 1 2 ) 9341,9342,9341 

9341 CY 12=(V2-V 1 >/b 12 
OX 1 2= < X 1 -X2 ) /D 1 2 

9342 I F < 023 >9343 , 9344 , 9343 

9343 0723= < V3-V2 > /D23 
0X23= < X2-X3 > /D23 

9344 IF <034 >9345, 9346, 9345 

9345 CV34= < V4-V3 >/D34 
0X34=<X3-X4>/034 

9346 I F < 04 1 >9347 . 9348 , 9347 

9347 CY4 1=<V 1-Y4 >/D4 1 
0X4 1 = <X4-X 1 >/D4 1 

0 05 COMPUTE MfiX LENGTH OF QUAD 

9348 ST=ftBS<X1-X3> 

ST2-SQ2F <X2,X4,Y2,V4, .0, ,0> 

ST=ftMftX 1 < ST , ST2 , 0 1 2 , 023 , 034 ,041) 

0 D. START LOOP OVER NULL PTS. 

342 KQ= 1 
L= 1 

343 l=KQ 

IF(IPS) 9360,9360,9350 
9350 1 F CL— 1 PS > 9355,9352,9352 
9352 IF(L- IPF> 9360,9360,9355 
9355 0 1 < JO >= . 0 
02<JC>=.0 
C3<JC>=.0 
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GU TO 541 
9360 I S= 1 

XCQ=T < I > 

VCQ=TCI+1> 

2CQ=T C I +2 > 

XNQ=T < I +3 > 

VNQ=T < I +4 ) 

344 ZNQ=T C I +5 > 

C E. COMPUTE DISTANCE BETWEEN QUAD AND NULL F'T . 

C . DETERM I N METHOD 

345 RFQ=S02F < XC , XCO , VC , VCO , ZC , ZCQ > 

I F < RPQ-ST*4 >350 , 350 ,460 

350 X=(XCO-XC >*XX+ C VCQ-VC >*YX+C2CQ-ZC >*2X 
S-CXCQ-XC >*XY+< VCO-VC >*VV+CZCQ-ZC >*ZV 
Z= CXCQ-XC >*XN+C VCQ-VC >*VN+ (ZCO-ZC >*2N 
IFCRPQ-ST+2 . 0 >355,355,400 

C F. COMPUTE UELGCITV COEF . BV EXACT METHOD 

355 R1=SQ2FCX , X 1 , V, V 1 . 2 , 0 . > 

R2=SQ2F < X , X2 , V , V2 , Z , 0 . > 

R3= SQ2F C X . X3 , V , Y3 , Z , 0 . > 

R4=SQ2F (X . X4, V , V4 . 2 , 0 . > 

IF CCFi1+F;2> .LE. D 12 > GO TO 1000 

IF <(R2+R3> .LE. 023) GO TO 1000 

IF CCR3+F;4> .LE. 034 > GO TO 1000 

IF <(F:4+Fi1 > .LE. D41) GO TO 1000 

CLA 1 =ALOO C <R 1 +R2-D 12 >/ < R 1 +R2+D 1 2 > > 

CLA2=AL0GC CR2+R3-D23 )/<R2+R3+023 > > 

CLR3=AL0G< CR3+R4-D34 >/CR3+R4+D34 > > 

CLR4=RL0G(<Fi4+R 1 — D4 1 >/<R4+Fl 1+D4 1 >> 

TUX=CV 1 2+CLA 1 +CY23*CLA2+CV34*CLR3+CV4 1*CLA4 
TUV=CX 12+CLfl 1 +CX23*CL A2+CX34 *CLR3+CX4 1*CLR4 
TUZ=0 . 

I FCABSCZ/ST )- . 0 10 > 375 , 36 1 , 36 1 

361 280=2**2 

E1=ZSQ+(X-X1>**2 
E2=2SQ+CX-X2 >**2 
E3=ZSQ+(X-X3>**2 
E4=ZSQ+(X-X4 )**2 
H 1= < V— V 1 >*CX-X1 > 

H2=(V-V2 >*<X-X2> 

H3=CY-Y3>*CX-X3> 

H4=CY-Y4 >*CX-X4 > 

IFCCI 12)363,363,364 

363 US 1 = < CM 1 2*E 1 -H 1 ) / C 2*F: 1 > 

US 2= < CM 1 2*E2-H2 > / C 2*Fi2 > 

AT 1=RTANCUS1 > 

AT2=ATAN<US2 > 

TUZ=AT 1-RT2 

364 IFCCI 23 >366 .366. 367 

366 AT3=RTAN < C CM23*E2-H2 > / C Z*R2 > > 

AT4=RT AN < < CM23*E3-H3 )/CZ*R3 > > 

T0Z=TU2+AT3-AT4 

367 IFCCI 34 >368 . 368 . 369 

368 RT5=HTANC CCI134*E3-H3 >/<Z*R3 > > 

RT6=ATRNC CCI134*E4-H4 >/CZ*R4 > > 

T0Z=T0Z+flT5-RT6 

369 IFCC 141 >370,370.375 

370 RT7=ATRHC COM4 1*E4-H4 >/<Z*R4 > > 

AT8=ATANC COM4 1*E 1-H 1 )/CZ*F; 1 > > 

TU2=TUZ+RT7-RT8 

375 GO TO 450 

C G. COMPUTE UELGCITV COEF. BV QUADRAPGLE METHOD 
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400 RPQ3=F:PQ**3 

RPQ7-<RPQ3**2)*RPQ 
US 1= X/RPQ3 
XSQ=X*+2 
VS0=Y**2 
ZSQ=Z**2 

PS=VSQ+ZSQ-4.*XSQ 
QS=XSQ+ZSQ-4 . ♦VSQ 
US2=X*<9. +PS+30 . *XSQ )/RPQ7 
US3=3 . *V*PS/RPQ? 

WS4=3 . *X*QS/RPQ7 

TUX=R*US 1-C I XV+WS3-C I XX*US2-C I VY+US4 
US1=V/RPQ3 

US2=V*<9 . *QS+30 . *YSQ)/RPQ7 
TW«fi*US 1-C I XX*HS3-C I XWIS4-C ! VV*US2 
TUZ=Z+ ( R /RPQ3-3 , + CCI XX*FS-5 . *C I XV*X*V+C I YV*QS ) /RPQ7 > 
450 UX( I S )=TUX*XX+TUY*XV+TUZ *XN 
UV< I S >=TUX+YX+TUV+YY+TU2*VN 
UZ( I S )=TUX*ZX+TW*ZV+TUZ*ZN 
GO TO 470 

C H. COMPUTE UELQC I TV COEF. BV MONOPOLE METHOD 

460 RRPQ3=fl/<RPQ**3> 

UX<lS>= < XCQ-XC )*nRPQ3 
UV ( I S >= < VCQ-VC >*RRPQ3 

UZ C I S )= < 2CQ-ZC X+HRPQ3 

C I . REFLECT HULL PT . IN PLRHE OF SVMETRV 

470 GO T0C480, 485, 490, 495, 500,505, 5 10, 5 15), IS 
C DO LOOPS SET UP TO FORCE USE OF INDEX REGISTERS 

480 J1=JU 
J2=JC 
UDX=UX < 1 ) 

UDY-UVC 1 ) 

UDZ=UZ( I ) 

U K J 1 )=UX( 1 > 

U2<J1 )=UX( 1 ) 

U3CJ1 )=UX( 1 > 

UKJ1+1 )=UV< 1 ) 

U2< J1+ 1 >=UV( 1 ) 

U3<J1+1 >=UV< 1 > 

U 1 ( J 1+2 >=UZ( 1 ) 

U2<J1+2)=UZC 1 ) 

U3<J1+2 >=U2< 1 ) 

IF(SVfl) 530,530,481 

481 IS=2 

C XZ SVMETRV 

VCQ=-VCQ 
GO TO 345 

485 I F <SVM— 1 >517,517,486 

C XV SVMETRV 

486 18=3 
ZCQ=-ZCQ 
GO TO 345 

490 ISM 
YCQ=-VCQ 
GO TO 345 

495 IFCSVM-2 >5 16, 5 16,496 

C VZ SVMETRV 

496 IS=5 
XCQ=-XCO 
GO TO 345 

500 18=6 
VCQ=-VCQ 
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GO TO 345 
505 IS=7 
2CQ=-2CQ 
GO TO 345 
510 IS=8 
VCQ=-VCQ 
GO TO 345 

J. ADD CONTRIBUTIONS OF RLL REFLECTIONS 

515 0 K0 1 >=U K J 1 >+UX<8 >+UX<? >+UX<6 >+UX<5 ) 

02 <0 1 )=02<0 1 )-UX<8 >+0X<7 >+0X<6 >-0X<5 ) 

03(0 1 >=03 <0 1 >-0X<8 >-0X<7 >+0X<6 >+0X<5> 

0 1 < 0 1 + 1 >=0 1< 0 1 + 1 >-0Y < 8 HUY < 7 HUY < 6 >-UY < 5 > 

02 < 0 1 + 1 >=U2 < J 1+ 1 >+UV < 8 HUY < 7 HUY < 5 HUY < 5 > 

U3( J 1+ 1 >=03<0 1+ 1 >+UV<8 HUY <7 HUY <6 HUY <5 > 

U 1 < J 1+2 >=U 1 < J 1+2 HUZ<8 >-UZ<7 >+U2<6 HU2<5 > 

U2< 0 1 +2 >=U2< 0 1 +2 HUZ < 8 )-UZ < 7 )+UZ < 6 >-UZ < 5 ) 

U3< J 1+2 >=U3< J 1+2 >+U2<8 )+U2<7 )+U2<6 >+UZ<5 ) 

516 U K J 1 >=U 1 < J 1 >+UX<4 >+UX<3 > 

U2CJ1 )=U2<J1 >+UX(4>-UX<3> 

U3CJ1 )=U3< J 1 )-UX<4 )-UX<3 ) 

U 1 < J 1 + 1 >=U 1 < J 1 + 1 )+W ( 4 )-UV ( 3 > 

U2< J 1+ 1 )=U2<J 1+ 1 >+UY<4 >+UV(3 > 

U3 <01+1 )=U3 <01+1 )-UV < 4 >+'0V < 3 > 

U 1 < J 1 +2 >=U 1 < J 1 +2 >-U2 < 4 >-U2 < 3 ) 

U2 <01+2 )=U2 < J 1 +2 >-U2 < 4 HU2 < 3 ) 

V3 < J 1 +2 >=U3 <01+2 )+UZ <4 )+UZ<3 > 

517 U 1 <0 1 )=U KOI )+UX<2 ) 

U2<J1)=U2<01 )-UX<2) 

03 < 0 1 )=U3 < 0 1 >+ UX < 2 ) 

U 1 < J 1+ 1 >=U 1 <01+1 )-UY<2 ) 

02<01+1 )=02<01+1 )+0Y<2> 

03 < 0 1 + 1 >=03 <01+1 >-0V < 2 > 

0 1 < 0 1 +2 >=U 1 < J 1 +2 >+U2 < 2 > 

02 <01+2 >=02 <01+2 >-0Z < 2 > 

03 < 0 1 +2 >=03 < 01+2 >+0Z < 2 ) 

530 C. 1 < J2 >=XNQ*0 KOI >+VMQ*U 1 <0 1 + 1 >+ZNQ*U 1 <0 1 +2 > 

C2<02 >=XN0*02<0 1 >+VM0*U2<0 1+ 1 >+ZMO*02<0 1+2 > 

C3C02 >=XNQ*03<0 1 >+YNQ*03< J 1+ 1 >+ZNQ+U3< J 1+2 > 

540 00=00+3 

541 OC=OC+ 1 

C D. WRITE COEFICIENTS 

C K. WRITE COEF. ON TAPE OR DRUM IF STORAGE AREA IS FULL 

545 I F < 00- 1 00 1 >570 ,555,555 
555 00=2 

OKI >=BLK 
02 < 1 >=BLK 
03 < 1 >=BLK 

I F<BLK-636 . 0 > 560,563,566 
560 WRITE <01> ELK, 01, 02, 03 
GO TO 56S 
563 REWIND 01 

566 MRITE< 1 1 > BLK,01,02,03 
568 BLK=BLK+ 1 . 

570 IF <00001 >580,571,571 

571 I DW= I DW+ 1 

WRITE <02 > I DU, Cl 
WRITE <08 > IDW.C2 
WRITE <03 > IDW,C3 
576 00=1 
580 KQ=KQ+7 
L=L+ 1 

C L. END OF LOOP OOEFl NULL PTS. 
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IFCKQ-k'M >343.343,581 
581 C 1 < JC >=0 
C2(JC)=0 
C3< JC )=0 

I F <KQ-KMM >54 1 , 585 , 585 

585 P=P+1 
K=K+7 
J=J+20 

I FCK-KM >586,586,600 
11. END OF LOOP ODER QUADS. 

M. READ NEXT BLOCK OF B RRRAV IF NEEDED 

586 I F< J-24 1 >296, 590 . 590 
590 READ <04 ><B< I ), 1 = 1,241 > 

J=2 

I F<B< 1 )-P >595, 296 . 595 
595 WRITE <6,98> B<1>,P 
STOP 

600 IF<BLK-636.0> 610,620,630 
C 0. WRITE REMAINING COEF. ON TAPE OR DRUM 

610 WR I TE<0 1 >BLK, U 1 , U2, 03 
REWIND 01 
GO TO 640 
620 REWIND 01 

630 WR I TE < 1 1 > BLK,U1,U2,U3 
REWIND 11 

640 WRITE <02 > IDW.C1 
WRITE <08> IDW‘C2 
WRITE <09 > IDW.C3 
REWIND 02 
REWIND 03 
REWIND 04 
REWIND 08 
REWIND 09 

C P. TRANSFER TO PFPS3 

GO TO 5000 

1000 WRITE<6,2000> L.P 

2000 FORMAT <3H L= ,15. 20X , 3H P= . F5 . 1> 

5000 CONT I NUE 
STOP 2 
END 

FUNCT I ON SQ2F (X 1 . X2 . V 1 . V2 . 2 1 , 22 > 

X=X 1-X2 
V=V 1-V2 
Z=21-Z2 

RS = Z !t,+ 2+V + *2+X* + 2 

R=ABS < X >+ ABS < V >+ABS < Z >+ 1 . 0E-20 

Fl=R+RS/R 

Fi= . 25+R+RS/R 

F;= R+RS/R 

SQ2F= . 25+R+RS/R 

RETURN 

END 
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C 

c 

c 

c 



PROGRAM PFP3(0UTPUT= 128, TRPE02, TRPE08, TRFE09, 

1 TfiPE 12, TRPE03, TRPE6=0UTPUT, TRPE2=TRPE02, 

2 TRPE8=TAPE08, TRPE9=TRPE09, TRPE3=TAPE03 ) 

XV2 POTENT I RL FLOW PROGRAM UERSION 1 SECTION 3 
SOLOES MATRIX EQUATION FOR SOURCE DENSITY 

COMMON SN < 654 ) , U I P ( 650 >, S ( 5 , 650 > , PROP ( 1 5 > , US <220 ) , DM < 650 ) , 

1 B < 220 > , COEF ( 900 ) , XN ( 650 ) , VN ( 650 ) , ZN ( 650 > 

EQU I URLENCE <US(213),EPS) , (KK,B(201)) 

EOU I URLENCE < M I X , US < 205 ) ) , CM I Y, US < 206 > ) , < M I Z , US < 20? ) ) 

EQU I URLENCE ( US ( 20 1 ).NP),( US < 208 > . I PS ) , < US < 2 1 2 ) , I PF ) . < US (211), KM 
KKM.MK) 

5 FORMAT (49H 1XVZ POTENTIAL FLOW PROGRAM SECTION 3, UERSION 4 ) 



URITE (6 


,5) 


READ 


(03 )(PROB( 1 >, 1 = 1.15) 


WRITE 


( 6, IQQIXPROBd ), 1 = 1, 11 


1001 FORMAT ( 1H0, 15R4 ) 


READ (03 


) CUSCI ). 1 = 1,220), 1 ED ITS. IEDIT4 


READ 


( 03 ) ( SK 1 P . SK I P , SK 1 P , XN ( 1 ) . VN ( 1 


D=- . 5/3. 


14159265 


RERD 


(03 )(DM( 1 ), 1 = 1 , NP ) 


K 1= 1 




K2=NP 




240 FORMAT ( 18H0CHRNGES IN PROS - 15R4) 


IF< IPS)1 


220, 1220, 1231 


-j 23 1 READ 


( 12 )< B(K),K=1 , 15) 


WRITE 


( 6, 240 )(B(K ), K= 1 , 15 ) 


READ 


( 12 )( B(K>,K= 1,220) 


RERD 


( 1 2 )SK 1 P 


READ 


( 12)SKIP 


K 1=1 PS 




K2= 1 PF 




R. SET CONDITIONS FOR FLOU OF -1 IN X 


1220 FX=- 1 




FV=0 




FZ=0 




NF=- 1 





C B. COMPUTE INITIAL APPROXIMATION TO THE SOURCE 

1240 DO 1250 K= 1 , NP 

U I P < K >=XN ( K )*FX+YN ( K )*FY+ZN < K)*FZ-DM(K) 

S(5,K)=-UIP<K>+. 1 1936 

C C. SET PARTIAL SUM UECTOR TO ZERO 

1250 SN(K)=0. 

SN(NP+1)=G. 

SN(NP+2)=0 . 

SN(NP+3)=U. 

SNCNP+4 )=0 . 

UR I TEC 6, 99?) FX,FV,FZ 
URITE (6,998) 

998 FORMAT < 27H0 1 TERAT I ON SUM OF CHANGES , 9X . 1 HA , 1 0X , 2HB 1 , 1 0X , 2HB2 ) 
IT= 1 



IC-5 

IF (IPS) 1260,1260,1255 

1255 RERD( 12) (S(5,K), K=1.KK) 

DO 1256 K=1,KK 

DO 1256 1=1,4 

1256 S( I ,K)=S(5,K) 

C D. START ITERATION 

1260 6RND=0 
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1261 

1262 



1263 

1264 
C 

C 

c 

c 



65 

67 

68 

69 

70 
C 

80 



C 

c 

1290 

C 

91 

92 

93 

94 



ICO 



110 

120 



IF CNF) 1261, 1262, 1263 

READ <02 > I DU, COEF 

GO TO 1264 

READ COS) I DU, COEF 

GO TO 1264 

READ <09)1 DU, COEF 

J=0 

D. READ FIRST BLOCK OF COEF 

E. START LOOP OUER QUADS. 

DO 1290 K= 1 , HP 

F. PICK UP SOURCE DENS I TV 
SF=SCIC,K> 

G. START LOOP OUER NULL PTS. 

DO 1290 KF- 1 , HP , 5 

I F< 0-900)80,65, 65 
IF GIF >67,68,69 
READ <02)1 DU. COEF 
GO TO 70 

READ <08)1 DU, COEF 
GO TO 70 

READ <09)1 DU, COEF 
J=0 

H. COMPUTE PARTIAL SUM FOR NEXT 5 PTS. 
SNCKP )=SN < KP )+COEF < J+ 1 )*SF 

SNCKP+ 1 )=SN<KF'+ 1 )+COEF < J+2 )*SP 
SNCKP+2 >=SH<KP+2 >+COEF< J+3 )*SP 
SNCKP+3 )=SN<KP+3 >+COEF< J+4 )*SP 
SNCKP +4 )=SNCKP+4 HC0EFCJ+5 )*SF 
U=.J+5 

J. END OF LOOP OVER NULL PTS. 

K. END OF LOOP OUER QUADS. 

CONT I HUE 

L. COMPUTE NEU SOURCE 
IF CNF) 91,92,93 
REWIND 02 

GO TO 94 
REWIND 08 
GO TO 94 
REWIND 09 
PASS=1 .0 
SUM=0 . 



DO 100 K=K 1 , K2 
SHCK)=< SN<K)+UIPCK) )*D 
TEST=ABSCSNCK )-S< I C , K ) ) 
SUH=SUI1+TEST 

IF CTEST .GT. EPS) PASS=-1.0 
CONTINUE 

IF CRASS .EO. 1.0) 

IF CIT.GE.MIX) 



IF < I ED ITS .EO. 
IT-IT+1 
I C= I C- 1 

IF < 1C .EQ. 0) 
DO 110 K=K 1 , K2 
SC IC:,K)=SNCK> 
SNCK)=0. 

GO TO 1260 
A=0 . 

B 1=0. 

B2=0. 



GO TO 180 
GO TO 180 
0) WRITE (6. 99) 



GO TO 120 



IT, SUM 



DR=0. 

D1=0. 
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02=0 . 

DO 140 K= K 1 , K2 
DS9=2*S( 1 , K >-SN(K )-S(2, K ) 

IF(DS9 .GT . 0. ) GO TO 122 
ft=fi+5(2,K>-S(1.K) 

DR=DR-DS9 
GO TO 125 

122 A =ft +S( 1,K)-S(2,K) 

DR=DR+DS9 

125 DS 1=S<4 . K >-S<3, K > 

D$2=S(3,K)-S<2,K) 

083=QS 1-DS2 
D8S=3(2, K )-S< 1 , K > 

DS5=DS2-DSS 
DS6=D8 1-DSS 
DS4=DS2-S<1,K)+SN<K) 

DS?=DS3*DS4-DS5*DS6 
DS8*DS6+DS5-DS4*DS3 
IF(DS7 .GT. 0. > GO TO 120 
B 1=B 1-DS 1+DS4 +002+006 
D1=D1-D£? 

GO TO 130 

128 B 1 =B 1 +03 1 *084-002*036 
01=01+057 

130 IF (OSS .GT. 0. ) GO TO 132 
B2=B2-D£ 1 *005+002*003 
D2=D2-DS8 
GO TO 140 

1 32 B2-B2+03 1 *035-002+003 
02=02+000 
140 CONTINUE 
fl=fl/DR 
E'1=B1/D1 
B2-B2/D2 

IF< IT EQ. 6) GO TO 155 

Afi=A-A3 

AA=fiBS(AR) 

IF (FIR .GT. .02) GO TO 148 
DO 145 K=K 1 . K2 

0(5, K )=R* < SN OO-S ( 1 , K ) )+S ( 1 , K ) 

145 S'NCK )=0 . 

WRITE (6, 6000) 

5000 FORMAT (29.X, 17HR EXTRAPOLATION ) 

GO TO 160 
14S BB1=B1-B31 
BE: 1=AB0(BB 1 > 

BB1=50.*BB1 

BB2=B2-BS2 

BB2=ABS(BB2) 

BB2=50 . +BB2 
BBB=AB8(B1) + AB0CB2) 

IF ( (BB1 .GT. BBB ) OR. (BB2 .GT. EBB) ) GO TO 1 
DO 150 K=K1.K2 

S ( 5 , K )=S < 2 , K )+B 1 * < 8(1, K )-0 ( 2 , K > )+B2+ ( SN ( K )-S ( 2 , K ) > 
150 8N(K)=0. 

WRITE (6, 7000) 

7000 FORMAT (2QX, 17HB EXTRAPOLATION ) 

GO TO 160 
155 DO 158 K=K 1 , K2 
0(5,K)=SN(K) 

15S 8H(K)=0. 

160 10=5 
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WRITEC&, 161 > ft, 61, 

161 FORMAT < 29X , 3E 1 2 . 3 > 
flS=fl 
BS 1=6 1 
BS2=E2 
GO TO 1260 

180 WRITE <6. 99) IT. SUM 
DO 182 K=K 1 , K2 
182 S(1.K>=SNCK> 

WR1TEC03) <S< 1, K>, K= 1 , HP > 

99 FORMAT <4X, 1 3, E 18. 5) 

997 FORMAT C13H0 X UELOCM TV= , F4 . 1, 15H V UELOC I TV=, F4 . 1 , 
1 15H 2 UELOC 1TV=,F4. 1 > 

IF(F2)1400, 1390, 1400 

C P 1 IF THIS WAS HOT LAST FLOW, SET FOR NEXT FLOW 

1390 F2=FY 
FV=FX 
FX=0 
M I X=M I V 

m;v=miz 

NF=NF+ 1 
GO TO 1240 
1400 REWIND 03 

C F'2 RERD IN PFPS4 RND TRANSFER TO IT 

STOP 3 
END 
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APPEND i X MJ ~ XVZPF SECTION PF4 



PROGRAM PFP4 < OUTPUT, TAPE6=0UTPUT, TAPE03, TAF'E0 1 , TAPE 1 1 , 
1 TAPE3=TRF'E03, TAPE 1=TAPE0 1 ) 



XVZ POTENTIAL FLOW PROGRAM UERSION 1 SECTION 4 
COMPUTES VELOCITIES AND PRESSURE COEFFICIENTS FOR 
POINTS ON THE BODV 



COMMON VX 1 < 650 ) , U V 1 < 650 ),VZ1< 650 > , VX2 < 650 ) , VV2 < 650 > , UZ2 <650 ) 

1 , UX3 < 650 ) , UV3 ( 650 ) , UZ3 < 650 ) , SI <650), S2<650), S3 <650) 

2 , X<650), V<650), Z<650), T4<650), T5<650), T6<650) 

3 , DM <650) 

D I MENS I ON PROS < 1 5 ) , US < 220 ) . CU 1 < 1 000 ) . CV2 < 1 000 ) , CU3 < 1 000 ) 
EQUIVALENCE (US <20 1 ), NF' ), <US<2 1 1 ),KM>, <US<217), VX'I ), <US<218),UVI ). 
1 <US<2 19 ) . UZ I ), < US <208 ), I PS ) . <US<2 12 ), I PF ) 

EQIJ I VALENCE <M I X , US<205 ) ), <M I V, US <206 )), <MIZ, US <207 ) ) 

5 FORMAT (49H0XVZ POTENTIAL FLOU PROGRAM SECTION 4, VERSION 4 ) 

URITE <6, 5) 

READ <Q3XPR0B<I ), 1=1, 15) 

100 FORMAT < 1H , 15R4) 

URITE <6, 100 XF'FlOBC I ) . I = 1 , 15 ) 

C A. READ PARAMETERS AND SOURCE 

READ <03 ) <US< I ), 1 = 1,220), IED1T3. IEDIT4 

READ <03 XX< I ), V< I ), 2< I ), T4< I ), T5< I ), T6< I ) , SK I P, 1 = 1, NF' ) 

READ <03 )<DM< I ), ! = 1,NP) 

D=-.5/3. 14155265 

READ (03XSK I ), 1 = 1 ,NP) 

READ <03 X SIX I ), 1 = 1 ,NP) 

READ < 03 ) < S3 < I ) , I = 1 , NP ) 

J=1 



K 1= 1 
K2=NP 

IF< I PS) 108, 108, 102 
102 K 1=1 PS 
K2= I PF 



103 BBFl= 1 . 



B . READ F 


1 RST BLOCK OF COEF. 




ITAF'E=01 






READ <01) BB 


,CV1 ,CV2,CV3 




IF <BE'-BBR) 


300, 120,300 




120 DO 125I=K1,K 


2 




VX 1 < 1 }=- 1 . 0 


-SKI >*T4< 1 ) 


/D 


VVKI )= 


-SKI >+T5< 1 ) 


/D 


VZKI )= 


-S 1 < 1 )+T6( 1 ) 


/D 


VX2< 1 )= 


-SIX' 1 >*T4< 1 ) 


/ D 


VVKI )=- 1.0 


-SK 1 )*T5< 1 ) 


/D 


VZ2< 1 )= 


-S2< 1 )*T6( ! ) 


/D 


UX3< 1 )= 


-S3< 1 )*T4< 1 ) 


/D 


VV3< 1 )= 


-S3< 1 )*T5< 1 ) 


/D 


125 VZ3< 1 )=- 1 . 0 


-S3< 1 )*T6< 1 ) 


/D 


C. SET UP 


LOOP OVER QUADS. 




JC=2 






D. PICK UP SOURCE 





130 S1J=S1<J) 
S2J=S2<J) 
S3J=S3< J ) 



C E. SET UP LOOP OVER NULL PTS. 

DO 180 JP=K 1 , K2 

C F. COMPUTE PARTIAL SUM FOFi 3 COMPONENTS OF 3 VELOCITIES 

VX 1 < UP )=VX 1 < JP )+S 1 J*CU 1 < JC ) 

VV 1 < JP )=VV 1 < JP )+S 1J+CV 1 < JC+ 1 ) 
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UZ 1 UP >=UZ 1 UP )+S 1J+CU 1 UC+2 ) 

0X2 < JP )=VX2(JP )+S2J*CU2(JC > 

W2< JP )=UV2< JP )+S2J*CU2(JC+ 1 > 

UZ2 < JP >=UZ2 < JP )+S2 J* CU2 < JC+2 > 

VX3 < JP )=VX3 < JP >+S3J+C'J3 ( JC ) 

UV3< JP >=VV3UP )+S3J*CU3< JC+ 1 ) 

UZ3 < JP >=UZ3 < JP )+S3J*CV3 < JC+2 ) 

JC=JC+3 

C G. REfiD MORE COEF. IF NEEDED. 

IF < JC- 1000 >180, 135 ,135 
135 JC=2 

IF(BBF;-635.Q> 140, 150, 155 
140 FiERD <01 > BB,CU1,CU2,CU3 
GO TO 160 
150 REWIND 01 

155 RERD (11) BE,CU1,CV2,CU3 
160 BBFl=BBR+ 1 . 

IF (BBR-BB) 300, 180,300 
C H. END OF LOOP OUER NULL PTS. 

ISO CONTINUE 
J= J+ 1 

C I. END OF LOOP OUER QUADS. 

I F< J-NP > 130 , 130,200 
200 I F ( BBFl-635 . 6 > 231, 23 1 , 232 

231 REWIND 01 
GO TO 233 

232 REWIND 11 

C K. EDIT THE UELOCITIES ETC. AND WRITE THEM ON TfiF'E 

233 WF; I TE <03 XVX 1(1 > , UV 1 < I ), UZ K I ) , I = 1 , NP ) 

WR I TE <03 XUX2< I ), UV2< I ), UZ2( I ) ,1 = 1, NP > 

WRITE <03X0X3(1 ), W3< I >,UZ3( 1 > , 1 = 1, NP) 

235 FORMRT < 1H 1 , 15R4, SH PAGE =,115) 

REWIND 03 



1 P=K 1+40 
I S=K 1 
I PRGE= 1 

IF < I ED I T 4 .EO. 1) 00 TO 293 
IF (MIX . LE . 0 ) GO TO 265 
242 FORMAT <8H0 X FLOW) 

240 FORMAT <4H FT . , 10.X, 1HX,9X, 1HV,9X, 1HZ, 13X.2HUX,8X,2H*JV,8X,2HUZ, 
1 5HABS.U, 7X, 2HCF' , 6X , 6HS0URCE , 4X , SHU NORMAL) 

245 FORMAT <1X, 13, 4X ,3F 10 .5. 4X ,3F 10 . 5, IF 13 . 5. 2F 1 1 . 5 , E 12 2) 

250 I F < I F-K2 >255 , 255 , 260 

C J. COMPUTE PRESSURE AND BBS. VALUE OF UELOC I TV 

255 WR I TE <6, 235 XPROB< I ), I = 1 , 15 ), I PAGE 

WRITE (6.242) 

WRITE (6,240) 

DO 257 1 = I S, IP 

USO-UX 1 < I )+*2+UV 1 < I )+*2+UZ 1 < I )+*2 

UM=(ABS(UX 1 < I ))+fiBS(UV1(l ))+fiBS(UZ1(l )>)+■. 79 

UH=v'M+VSQ/UM 

UM= . 25+UM+VSQ/UM 

UM =.5*(UM+US0/UM) 

CP = 1 . -VSQ 

UNF: =UX 1 ( I )*T4 ( I ) +VV 1 < I )*T5 ( I ) +VZ 1 < I )+T6 ( I ) 

257 WRITE (6,245) I , X( I ), V< I ) , Z< ! ), UX 1 ( I ), W 1(1 ), UZ 1 < 1 ), VM 
1 , CP ,S1( I >,UNF; 

I S= I S+-50 
| p= l P+50 
I PROE= I PRGE+ 1 
IF(K2-IS) 265,260,250 
260 IP=K2 



10X, 
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GO TO 255 
265 IP=K1+49 
I S=K 1 

IF (HIV LE 0 ) GO TO 280 
26? I F ( I P-K2 >275 , 275 , 270 
270 I P=K2 

275 WR I TE <6, 235 ><PROB< I > , I = 1 , 15 > , I PAGE 

WRITE (6,27?) 

27? FORMAT (SHO V FLOW) 

WRITE (6,240) 

DO 278 1= IS, IP 

080=0X2 < I )++2+UV2( I )*+2 +022 < I >**2 

UM=(ABS(UX2( I > )+ABS(UV2( I > >+ABS(UZ2( I ) ) >* . 79 

UM=UM+USQ/UM 

UM= . 25*UM+USQ/UM 

OH = .5*(UM+VSQ/Uri) 

UNR =0X2 < I )*T4 ( I > +UV2 < I )*T5( I > +UZ2 ( I )*T6 ( I > 

CP = 1 . -OSO 

278 WRITE (6,245) I ,X( I ),Y( I ),Z( I ),UX2( I ),UV2( I >,U22( I ),0M 
1 .CP .82(1), OMR 
18=18+50 
IF- 1 P+50 
IPRGE-IPAGE+1 
I F( 18-1(2)267,267,280 
280 IF-K1+49 
I S=K 1 

IF (MIZ . LE . 0 ) GO TO 293 
282 I F ( I P-K2 >290 , 290 , 285 
285 IF-K2 

290 WRITE ( 6 , 235 ) ( P ROB ( I ) . I = 1 , 1 5 ) , I PAGE 

WRITE (6.292) 

292 FORMAT (SHO Z FLOW) 

WRITE (5,240) 

DO 291 1=18, IP 

US0=UX3( I )++2+UV3( I )**2+UZ3( I >+*2 
UM=(RBS(UX3( I ) )+ABS(UV3( I ) )+ABS(UZ3( I ) > )+ . 79 
UM=UM+USQ/UM. 

UM= . 25+UM+USQ/UM 
UM = . 5* ( UM+USQ /WM ) 

CP = 1 . -OSO 

UNR =UX3 ( I )*T4 ( I ) +UV3 ( I ) :t: T5 ( I ) +UZ2 ( I X+T6 ( I ) 

291 WRITE (6,245) I ,X( I ),Y( I ),Z( I ),UX3< I >,UY3< I ),UZ3( I ),UM 
1 ,CP ,83(1), UNR 

1 8= 1 8+50 
I P= I P+50 
I FfiGE= I F'RGE+ 1 
I F< I S-K2 )282 , 282 . 293 
C L. CHECK FOR fl FOURTH FLOW 

293 J = 1 

294 UX.I = WS(.J) 

UVi = WS(J+20> 

UZI = W8(J+40) 

295 IF (UX I ++2+UV I **2+02 1 *+2 ) 400,400,301 
301 I S=K 1 



IP=K1+49 

320 IF (IP — 

C N. 




325 I P=K2 
330 WRITE 



(6,235)(PR0B(I ), 1 = 1, 15), I PAGE 

(6,315)0X1, UVI, UZI 

(6,340) 



WRITE 

WRITE 



C 



n. COMPUTE UELOC I TV AND PRESSURE FuFi FOURTH FLOW 



DO 333 l = IS, IP 

UX4 =- (MX I *UX 1 < I HUY I *UX2 ( I HUZ 1 *UX3 ( I > > 

UV4 =- < UX I *UV 1 ( I )+UV I *UY2 ( I )+UZ I *U V3 ( I ) ) 

UZ4 =- (MX I *UZ 1 ( I )+UV I *UZ2 ( I >+UZ I *UZ3 (I >) 

MSQ=MX4 ++2+UY4 ++2+UZ4 **2 

Mf1=(flBS(MX4 >+RBSCUY4 )+ABS(UZ4 ) )* . 70 

UM=UM+USQ/UM 

MM= . 25+uri+MSQ/uri 

UM =.5*(UM+USQ/UM) 

CP = 1 . -(USQ >/( USU++2+UZI**2+UXI**2) 

333 WRITE (6,345) I , X ( I ) , V ( I ) , Z ( I ) , UX4 , UV4 , MZ4 , UM 
1 ,CP 
I $= 1 3+50 
I P= I P+50 
I PflGE= I PBGE+ 1 



I F(K2- 1 S >350, 

350 J = J+1 
GO TO 294 

3 1 5 FORMAT ( 1 9H0 ONSET FLOW , MX I = , F6 . 3 , 2X , 4HUV I = , F6 . 3 , 2X , 4HMZ I ■ , F6 . 3 ) 
340 FORMAT (4H FT., 10X, 1HX,QX, 1HV,QX, 1HZ, 13X, 2HUX , 8X , 2HUV , SX , 2HMZ , 1 OX , 
1 5HABS.U, ?Xj 2HCP ) 

345 FORMAT ( IX. 13, 4X, 3F 10 . 5, 4X,3F 10 .5, IF 13 . 5, IF 11 . 5 ) 

400 CONTINUE 



GO TO 5000 
300 CONTINUE 

WRITE (6,310) ITAPE,BBR,BB 

310 FORMAT (6H1TRPE , 12, 17H OUT OF F'OSS IT ION/ 1 14, 6F6 . 1 ) 
5000 CONTINUE 



STOP 4 
END 
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RPPEND I X U - XVZPF SECT l OH PF5 



C 

c 

c 

c 

c 



PROGRAM PFP5C I NPUT= 128, OUTPUT® 128, TAPE03, TAPE04, 

1 TAPE5= 1 NF'UT, TAPE6=0UTPUT, TRPE3=TAPEQ3 , TAPE4=TAPE04 ) 

XVZ POTENTIAL FLOW PROGRAM UERSION 4 SECTION 5 
COMPUTES UELOCITIES AND PRESSURE COEFFICIENTS FOR 
OFF BODY POINTS 

COMMON B<24 1 >, XP(500 >, YPC500 > , ZPC500 >, UX 1 <500 ), US <220 ) 

1 , UY 1 <500 ), UZ 1 <500 ), UX2<500 ), UY2<500 ), UZ2<500 ), UX3<500 ), W3<500 ) 

2, UZ3 < 500 ) , UX < 8 ) , UY < 8 ) , UZ<8 ).S 1 <650 ) , S2<650 ), S3 <650 ). U 1 <3 ), U2<3 >, U3 
3<3 > , PROB< 15) , XN<650), VN<650 >, ZN<650 >, TA<550 >, TX<650 ) 
4,TV<65u),TZ<650) 

EQU I UALENCE < KM . US <2 1 1 ) ) , < KM , MK ) , < NP , US < 20 1 ) ) , < SVM , US < 2 1 0 > ) 

1 , <V2, Y3 ) , <M I X, US <205 ) ), <M I V, US<206 ) ), <M I Z , US <207 ) ) 

INTEGER PAGE 
INTEGER SVM 

5 FORMAT <49H0XVZ POTENTIAL FLOU PROGRAM SECTION 5, UERSION 4 ) 

UP.ITE <6, 5) 

A. READ INPUT 

READ <5, 25) NOBP, I EDITS, I READ 

A. READ THE OFF BODY POINTS 
NOB=NOBP 

DO 10 1=1. NOB 

READ <5, 26) XP< I ). VP< I ),ZP< I ) 

IF <E0F<5) ,EQ. 0.) GO TO 10 
NOBF- 1 - 1 

URITE<6,9) NOBP. MOB 

9 FORMAT <1 HO, 1 5 , 3 1 H OFF BODY POINTS SPECIFIED NOT ,13) 

GO TO 11 

10 CONTINUE 

11 CONTINUE 

25 FORMAT <3 1 4 > 

26 FORMAT <3F 12.5) 

P=1. 

READ <03 > <PF;OB< I ). 1 = 1, 15) 

UR ! TE <6 . 90 > <F'ROB< I > . I = 1 , 15 > 

90 FORMAT < 1H0, 18A4) 

UR I TE<6, 9 1 ) NOBP , I ED I T5 , I READ 

91 FORMAT <8H0N0E'P =, 14 /8h’ I ED I T5= . 14/SH I READ =, 14) 

UR I TEC 6, 92) 

92 FORMAT < 17H0 OFF BODY POINTS / 4H FT., 1 1X, 1HX, 12X, 1HV, 12X, 1HZ) 

UR I TE<6 . 93 ) < I ,XP< I ),VF'< I ),ZP< I ), 1 = 1, NOBP ) 

93 FORMAT < IX, 1 3, 2X , 3F 13.5) 

B. READ THE PARAMETERS. T ARRAY AND SOURCE FROM TAPE 31 

READ <03) < US < I ), 1 = 1,220) 

READ <03 ) <TX< ! >, TV< I ), TZ< ! ), XN< I ), VN< I ), ZM< I >,TA< I ), 1 = 1, NP ) 

READ <03 ) SKIP 

FORNERLV : US <220) .EQ. 2. 

I F < US < 220 ) . EQ . 5 . ) READ < 03 ) SK I P 
READ <03 ) <S1< I ), l = 1,NF ) 

READ <03 > <S2< I ), 1 = 1, NP > 

READ <03 ) CS3< I ), l = 1,NP) 

0. READ THE FIRST BLOCK OF THE B RFlRAV 
READ <04 ) <B< I ), 1 = 1,241) 

K= 1 
J=1 

DO 100 1=1, NOBP 

D. SET THE PARTIAL UELOCITY TO THE FREE STREAM UELOCITV 
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UX1(I)=-1.0 
WKI >=0. 

UZ K I )=0 . 

UX2< I )=0. 

UY2CI >=-1.0 
UX3( I >=0 . 0 
W3<l >=0.0 
UZ3< I )=-1 .0 
100 UZ2< I )=0. 

290 I F 1 >-P >29 1 , 295, 29 1 

291 WRITE <6,98) B< 1 >,P 

98 FORMAT <28H0 POINTS OUT OF ORDER B< 1 >=. 1F4 0,4H P=,1F4.0> 

STOP 

C E. STRRT LOOP ODER THE QUfiDS. 

295 J=2 

C FI PICK UP QURD. INFORMATION 

296 X1=B(J> 

V 1=B<U+ 1 > 

X2=B<J+2> 

V2=B< J+3) 

X3=B< J+4 > 

X4=B<J+5> 

V4=B< J+6 > 

XC=TX(K > 

VC=TV<K> 

ZC=TZ<K> 
ft =Tfi(K> 

XX=B(J+7> 

VX=B<J+S> 

ZX=B< J+9 > 

XV=B< J+ 10 > 

VV=B< J+ 1 1 > 

ZV=B<J+ 12) 

C F2 COtIPUTE LENGTH OF SIDES OF OURD. 

D 12=SQ2F<X 1 , X2, V 1 , V2, 0. , 0 . > 

D23=SQ2F <X2,X3,V2,V3, .0, .0) 

D34=SQ2F<X3, X4 . V3 , V4, . 0, . 0 > 

D4 1=SQ2F<X4 .XI. V4 - . V 1 . .0. .0 > 

C F3 COMPUTE SLOPE OF SIDES 

I F < X2-X3 >305 , 300 .305 
300 Cl 23=1. 

GO TO 310 

305 Cf123= < V2- V3 ) / < X2-X3 > 

Cl 23=0. 

3 10 I F < X3-X4 >3 15,311,315 

311 Cl 34=1. 

GO TO 320 

3 15 CM34=(V4-V3 >/<X4-X3 > 

Cl 34=0. 

320 I F<X4-X 1 >325, 32 1 , 325 

321 Cl 4 1=1. 

GO TO 330 

325 CM 4 1=<V 1-V4 >/<X 1-X4 > 

C 14 1=0. 

330 I F <X 1-X2 >335 , 33 1 . 335 

331 Cl 12=1. 

GO TO 340 

335 CM 1 2= < V2-Y 1 >/<X2-X1 > 

Cl 12=0. 

C F4 COMPUTE OURDRRPOLE MOMENTS 

340 CIXX=B<J+1?> 

CIXV=B(J+18> 
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C1W=B(J+19> 

F5 COMPUTE SIN FIND COS OF SLOPE ANGLE FOR EACH SIDE 
CV12=(V2-V1)/D12 
CV23=(V3-V2>/D23 
CY34= < V4- V3 > /D34 
CV4 1=(V 1-V4 >/D4 1 
CX 1 2= ( X 1 -X2 > /D 1 2 
CX23=(X2-X3>/D23 
CX34=(X3-X4>/D34 
CX4 1=<X4-X 1 >/D4 1 
C F6 COMPUTE MAX D I AG INAL 

ST=SQ2F(X 1 , X3, V 1 , VS , 0 . , 0 . ) 

ST2=SQ2F(X2 . X4, V2, V4, .0. . 0 ) 

I F (ST-ST2 >34 1 , 342 . 342 

341 ST-ST2 

C G. START LOOP DOER THE OFF BODV POINTS 

342 DO 530 JQ=1,N0BP 

IS-1 

XCC|=XP(JQ> 

VCQ=VP( JO ) 

2CQ=2F'( JO > 

J1=1 

345 RPQ-SQ2F ( XC , XCQ , VC , VCQ , 2C , ZCQ > 

C H. DETERM I N METHOD 

I F(RPQ-ST*4 >350 , 350, 460 
350 X= (XCQ-XC >*XX+<VCQ-VC >*VX+(2C0-2C >+2X 
V= (XCQ-XC >*XV+(VCQ-VC >*W+(ZCQ-ZC >*2V 
2= ( XCQ-XC >*XN < K >+ ( VCQ- VC >* VN <K )+ < ZCQ-2C >*2N < K > 
IF(RPQ-ST*2 . 5 >355, 355 , 400 

C I. COMPUTE I MOUSED UELOC I TV BV EXACT METHOD 

355 R 1=SQ2F<X , X 1 , V, V 1 , 2 , 0 . > 

R2=SQ2F«’X2.V,V2,z',0. > 

R3-SQ2F ( X , X3 , V , V3 ,2,0. > 

R4=SQ2F ( X . X4 , V , V4 , 2 , 0 . > 

IF((R1+R2>.LE.012> GO TO 1000 
I F ( <R3+R2 > . LE . D23 > GO TO 1000 
I F ( ( R3+R4 > . LE . D34 > GO TO 1000 
IF((R1+R4>.LE.D41> GO TO 1000 
CLfl 1-RL0G«R 1+R2-D 12 >/<R 1+R2+D 12 > > 

CLR2=AL0G( (R2+R3-D23 >/(R2+R3+D23 > > 

CLR3=RL0G'( (R3+R4-D34 >/(R3+R4+D34 > > 

CLR4-AL0G< CR4+P, 1-D4 1 >/<R4+R 1+D4 1 > > 

TUX=CV 12*CLA 1+CV23+CLA2+CV34+CLA3+CV4 1*CLR4 
TUV=CX 12+CLA 1+CX23*CLA2+CX34*CLA3+CX4 1+CLR4 
T02=0. 

I F(ABS<Z >- . 00 1*ST >375,35 1 . 35 1 
361 230=2**2 

E1=2SQ+(X-X1 >+*2 
E2=2SQ+ ( X-X2 >**2 
E3=2SQ+ CX-X3 >**2 
E4=2SQ+<X-X4>**2 
H 1 = ( V- V 1 )* ( X- X 1 > 

H2= ( V-V2 >* ( X-X2 > 

H3= < V-V3 >* (X-X3 > 

H4=(V-V4 >+<X-X4 > 

IF (Cl 12 >363, 353. 364 

363 US1=(CM 12+E 1-H 1 >/<2*R 1 > 

I4S2= ( CM 1 2+E2-H2 >/(2*R2 > 

AT 1= AT AN (MSI > 

AT2=RTRN(US2 > 

TV2=AT 1-RT2 

364 IF (Cl 23 >366, 366, 367 



366 fiT3=RTfiN< CCN23+E2-H2 >/<Z*R2 > > 

RT4=flTRN< <CM23*E3-H3 >/<Z*R3 > ) 

TUZ=TUZ+RT3-RT 4 
357 I F<C 1 34 >363,368,369 

368 RT5=RTflf« <Ct134*E3-H3 >/<Z*R3 > > 

RT6=RTRN< <CM34*E4-H4 >/<Z*R4 > > 

TUZ=TUZ+RT5-fiT 6 

369 IF(C 14 1)370,370,375 

370 RT7=RTRNC(CM4 1+E4-H4 >/CZ*R4 > > 

RT8=RTftH< COM4 1*E 1-H 1 >/<2+R 1 > > 

TUZ=TUZ+fiT7-fiT8 

375 GO TO 450 

C J. COMPUTE I MOUSED UELOC I TV- BV QURDRRPOLE METHOD 

400 RPQ3=RPQ**3 

RPQXRPQ3**2 >*RPQ 
US 1= X/RPQ3 
XS0=X**2 
V S0=V**2 
ZSQ=Z**2 

PS=VS0+ZS0-4 . +XSQ 
OS=XSQ+ZSO-4 . *VSQ 
WS2=X*<9 . *PS+30 . +XSQ >/RPQ? 

1-133=3 . +V+PS/RPQ7 
US4=3 . *X*QS/RPQ? 

TUX=fl*US 1-C I XV+US3-C I XX*US2-C I VV+US4 
US 1=V/RPQ3 

US2=V*<9 . *08+30 . *VSQ J/RP07 
TUV=fi x ;,6 1-Ci XX+US3-C ! XV*US4-C i VV+US2 
TUZ=Z* < ft /RP03-3 . *<C I XX+PS-5 . *C I XV+X+V+C I VV+OS J/RPQ7 > 
450 UX< I S >=TUX*XX+ T U V*XV+ TUZ+XH < K > 

W< I S >=TUX*VX+TUV*VV+TUZ*VH<K > 

02 C ! S >=TUX*ZX+ TU V+ZV+TUZ+2M ( K > 

GO TO 470 

C K. COMPUTE I HOUSED UELOC I TV BV MOHOPOLE METHOD 

450 RRPQ3=R/(RPQ**3> 

UX < I S >= < XCQ-XC >+RPiP03 

UV(IS>= (VCO-VC >*RRP03 
UZ(IS>= CZCQ-ZC >+RRFQ3 

C L. REFLECT OFF BODV POINT IN PLRME OF SVMETRV 

470 GO T0C480, 485, 490. 495, 500, 505, 5 10. 5 15), IS 

480 U 1 (J 1 )=UX< 1 > 

U2<ji >=uxx i > 

U3(.Ji>=uxa:> 

UKJ1+ 1 >=UV( 1 ) 

U2< J1+ 1 )=UY< 1 > 

U3(J1 + 1 )=UV< 1 > 

U1<J1+2>=U2< 1 > 

U2(J1+2)=U2<1> 

U3(J1+2)=UZ( 1 ) 

IF(SVM) 525,525,481 

481 IS=2 

C XZ SVMETRV 

veq=-vco 

GO TO 345 

485 i F < SVM— 1 >517.517,486 

C XV SVMETRV 

486 I S=3 
2C0=-ZCQ 
GO TO 345 

490 I S=4 

VCQ=-VC:0 
GO TO 345 
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495 1 FCSVM-2 >5 16, 5 16, 40fi 

C VZ SVMETRV 

496 I S=5 
XCQ=-XCQ 
GO TO 345 

500 IS=6 
VCQ=-VCQ 
GO TO 345 
505 IS- 7 

ZCQ=-ZCQ 
GO TO 345 
510 IS=8 

VCQ=-VCQ 
GO TO 345 

C M. ADO CONTRi BUT IONS OF ALL REFLECTIONS 

5 15 0 K J 1 >=0 K J 1 >+0X(8 >+0X(7 >+0X(6 >+0X(5 ) 

02(0 1 )=02(0 1 >-0X(8 >+0X(7 >+UXC6 )-UX(5 > 

03 ( J 1 >=03 < J 1 >-0X ( 8 >-UX ( 7 >+0X ( 6 >+0X ( 5 > 

0 1 < J 1+ 1 >=U 1< J 1+ 1 >-0V(8 >+0V < 7 >+U V ( 5 >-UYC5> 

02 <01+1 >=02 (01+1 >+0V < 8 >+0V < 7 >+UV ( 6 >+0V ( 5 > 

03 < J 1 + 1 >=03 < J 1 + 1 HOY < 8 >-0V ( 7 >+0V ( 6 >-0V < 5 > 

0 1 ( J 1+2 >=0 1 (01+2 >-0Z(8 >-0Z<7 >+UZ(6 >+0Z(5 > 
02(01+2 >=02 (J 1+2 >+0Z < 8 >-0Z < 7 >+0Z ( 6 >-0Z(5 > 
03 ( J 1 +2 >=03 (01+2 >+0Z ( 8 >+0Z ( 7 >+0Z ( & >+0Z ( 5 > 

5 16 0 1 (0 1 >=0 1(01 >+0X(4 >+0X(3 > 

02(0 1 >=02(0 1 >+0X<4 >-0X<3 > 

03(0 1 >=03(0 1 >-0X(4 >-0X<3 > 

0 1 < 0 1 + 1 >=0 1 ( 0 1 + 1 >+UV ( 4 >- 0 V ( 3 > 

02 (01+1 >=02 (01+1 >+0V ( 4 >+0V ( 3 > 

03(01+ 1 >=03(01+ 1 HUY <4 >+0V<3 > 

01(01+2 >=0 1(01+2 >-0Z ( 4 >-0Z ( 3 > 

02(01+2 >=02 ( 0 1 +2 >-UZ ( 4 >+0Z ( 3 > 

03(0 1 +2 >=03 ( J 1 +2 >+0Z ( 4 >+0Z ( 3 > 

517 01(01 >=01(01 >+0X(2> 

02(01 >=02(01 >-0X(2> 

03(01 >=03(01 >+0X(2> 

0 1(01+1 >=0 1 (0 1+ 1 >-UV(2 > 

02(01+1 >=02(01+1 >+0Y(2> 

03 < 0 1 + 1 >=03 < 0 1 + 1 >-0Y ( 2 > 

01 ( 01+2 >= 0 1 ( 01+2 >+ 0Z ( 2 > 

02(01+2 >=02 (01+2 >-0Z ( 2 > 

03 ( 0 1 +2 >=03 ( 0 1 +2 >+0Z ( 2 > 

525 L=P 

OX 1 (OQ >=0X 1 (00 >+01(1 >*£ 1 (L > 

OY 1 ( OQ >=0V 1 < 00 >+0 1 ( 2 >+S 1 < L > 

OZ 1 (OQ >=0Z 1 (OQ >+0 1 (3 >*S 1 (L > 

0X2 (OQ >=0X2 (OQ >+02 ( 1 >*S2(L > 

0V2 ( 00 >=0V2 ( 00 >+02 ( 2 >+S2 ( L > 

0Z2 ( OQ >=0Z2 < OQ >+02 ( 3 >*S2 ( L > 

0X3(00 >=0X3(00 >+03 ( 1 >*S3(L > 

0Y3 ( OQ >=0V3 ( 00 >+03 ( 2 >+S3 ( L > 

023 ( OQ >=023 < 00 >+03 ( 3 >+S3(L> 

530 CONTINUE 

C N. END OF LOOP GOER OFF BODY POINTS 

585 P=F'+ 1 
K=K+ 1 
0=0+20 

IFCK-NP >586,586,599 

586 IF (0-24 1 >296, 590,590 
0. REfiD NEXT BLOCK OF B ARRAY IF NEEDED 

590 READ (04 > (B( I >, 1=1,241 > 

0=2 



C 



c 



IF(B( 1 )-P)29 1,296, 291 
P. END OF LOOP O'JER QUADS 
599 CONTINUE 
PAGE = 1 

IF ( I EDITS .EO. 1) GO TO 825 

601 FORMAT (4H PT. , 1 IX, 1HX, 12X, 1HV, 12X. 1HZ. 14X.2HUX. 1 1X,2HUY, 1 1X,2HU 
1, 14X,2HCP) 

602 FORMAT (7H X FLOW) 

603 FORMAT (7H V FLOW) 

604 FORMAT (7H 2 FLOW) 

605 FORMAT (1 HI, 15R4, 10X, 15H0FF BODV POINTS , 10X,5HF'AGE ,13) 

IF (MIX.EQ.O) GO TO 700 
WRITE (6. 605) PFlOB , PAGE 
WRITE (6,602) 

WRITE (6,601) 

L I NE= 1 
LAST=53 

606 I F ( NOBP-LAST )607 ,610,610 

607 LAST =NOBF’ 

610 DO 615 l=LINE,LAST 

6 1 1 FORMAT < 1 X , 1 1 3 , 2X , 3F 1 3 . 5 , 2X , 3F 1 3 . 5 , 3X , F 1 3 . 5 ) 

C Q. COMPUTE PRESSURE AND EDIT 3 BASIC FLOWS 

CP 1= 1 . -<UX 1 ( I )**2+UV 1 < I )**2+U2 1 < I )**2 ) 

615 WRITE (6.611) I , XP( I ), VP( I ),ZP( I ),UX 1(1), UV1 ( I ), U21( I ), CP 1 
LINE=LAST+ 1 
LAST=L IHE+54 
PAGE=PAGE+ 1 

I F(L I NE-NOBF' >620 , 620,700 

620 WRITE (6, 605) PROB.PRGE 
WRITE (6,601) 

GO TO 606 

700 IF (MIY.EQ.O) GO TO SOO 
WRITE (6, 603) PROB.PRGE 
WRITE (6,603) 

WRITE (6,601) 

L I NE= 1 
LAST=55 

706 I F ( NOBP-LAST )707 . 7 1 0 . 7 1 0 

707 LRST=NOBP 

710 DO 715 INLINE, LAST 

CF'2=1 ,-(UX2( I )*+2+UV2( I )*+2+U22( I )**2 ) 

715 WRITE (6,611) I , XF’C I ), YP< I >.2P( I >, UX2C I ), UY2( I ),UZ2< I ),CP2 
LINE=LRSt+1 
LRST=L IHE+54 
PROE=PRC€+ 1 

I F(L I NE-NOBF' >720,720,800 

720 WRITE <6, 605) FRGB , PAGE 
WRITE (6,601) 

GO TO 706 

800 IF (MIZ.EO.O) GO TO 825 
WRITE (6, 605) PROB.PRGE 
WRITE (6,604) 

WRITE (6,601) 

L I ME= 1 
LRST=55 

806 I FCNOBP-LRST )807, 8 10, 810 

807 LAST=NOBP 

810 DO 815 I -LINE, LAST 

CP3= 1 . -<UX3( I >**2+UY3( I )**2+UZ3( I >**2 > 

815 WRITE (6,611) I , XP( I ) , VP( I ), ZP( I ),UX3( I ) , UV3( I ), UZ3( I ) , CPS 
LINE=LRST+1 
LAST=L IHE+54 



PRGE=PRG'E+ 1 

I F<L I NE-NOBP >820.. 820, 825 
820 WRITE (6,605) PROS, PAGE 
WRITE (6,601) 

00 TO 806 

825 J = 1 

826 IF ( I READ . EQ . 0 ) GO TO 827 
RERD(5,26> UX4,UV4,UZ4 

IF (E0FC5).NE. 0.) GO TO QOO 
00 TO 828 

827 UX4=US(J) 

I.JV4 = WS<J+20) 

UZ4 = US(J+40) 

828 CP<>X4**2+UY4*+2+UZ4+*2 
I F(CP >900,900, 830 

C R. COHPUTE FOURTH FLOW RND EDIT IT 

830 L I NE= 1 
LRST=5 1 

WRITE (6, 605) PROS, PAGE 

831 FORMAT C19H00NSET FLOW UX =,F7.3/15X,4HW =, 
1F7.3/15X.4HUZ =,F7.3) 

WRITE (6,831) UX4,UV4,UZ4 
WRITE (6,601) 

835 I F ( NQBP-LAST )83? , 840 , 840 

837 LAST=NGBP 

840 DO 845 I =L I HE, LAST 

UXF--UX4+UX 1 ( I >-UY4+UX2( I )-UZ4*UX3( I > 

UVF --UX4+UV 1 ( I )-UV4+W2< I >-UZ4+UV3< I ) 

VZF--VX4*VZ1< I >-W4*VZ2( I >-UZ4+VZ3< I ) 

CP4= 1 . - ( UXP**2+ WP**2+UZP+*2 ) /CP 
845 WR I TE (6,611) I , XP ( I ) , VP ( I ) , ZP ( I ) , UXP , UVP , UZP , CP4 
LINE=LRST+1 
LAST=LINE+54 
PAOE=PAOE+ 1 

I F ( L I NE-NOBP )850 , 850 , 860 
850 WRITE(6,605) PROB, PAGE 
WRITE (6,601) 

GO TO 835 
860 J = J+ 1 
00 TO 826 

1000 WR I TEC 6 .1001) JQ , L , XPC JQ > , VP( JQ ) . ZPC JQ ) 

1001 FORMAT (16H00FF BODY POINT , I3,23H OH BOUNDARY OF OURD 

1 3H X=, F 12. 5,5X, 2HV=, F 12.5, 5X, 2HZ=,F 12.5 > 

GO TO 530 
900 CONTINUE 

C S. REWIND TAPES AMD STOP 

REWIND 03 
REWIND 04 
STOP 5 
END 

FUNCT I OH SQ2F (XI, X2 , V 1 , V2 , Z 1 , Z2 ) 

X=X1-X2 
V=V 1-V2 
Z=Z1-Z2 

RS=Z**2+Y**2+X**2 

R=ABS(X )+RBS(V)+ABS(Z )+ 1 .0E-20 

Fi=Fi+RS/R 

R=.25*R+R8/R 

Fl= R+RS/R 

SQ2F= ,25*R+RS/R 

RETURN 

END 
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APPENDIX VI - XVZPF SECT 1 OH PF6 



PROGRAM FFF’6 ( I NPUT= 1 28 , TAPE 1 6 , 0UTPUT= 1 28 , TAFE03 , TAPE04 , 

1 TAPE5= I NPUT TAPE6=0UTPUT , TRPE3=TRPE03 , TftPE4=TAPE04 ) 

XVZ POTENTIAL FLOW PROGRAM VERSION 4 SECTION 6 
COMPUTES VELOCITIES AND PRESSURE COEFFICIENTS FOR 
OFF BODV STREAMLINES 



C 



C 



C 



VERSION 4 > 

, 14 , 1GH STEPS OF 



COMMON XP < 1 00 >, VP ( 1 00 ) , ZP <1 00 ) , VX 1 < 1 00 ) , WS ( 220 > 

1 , VV K 1 00 ) , VZ 1 < 1 00 ) , VX2 <1 00 > , V V2 <100), VZ2 < 1 00 ) , VX3 < 1 00 ) , VV3 < 1 00 ) 

2 . VZ3< 100 ), VX<8 ) . UV<8 ), VZ<8 ) . S 1 <650 ) . S2<650 ), S3<650 ) , V 1 <3 ) . U2<3 ) , U3 
3<3) , PROB < 1 5 ) , XN < 630 ) , VN < 650 ) , ZN < 650 ) , TX < 650 ) , TV < 650 > , TZ < 650 ) 

1 , B < 1 3000 > , ' XT < 100 ), VT< 100 ), ZT< 100 ) . AP<5 >, GM < 4 > , SKY < 1 00 ) , 

1 SKZ < 1 GO > , SKX < 1 00 ) , DX < 1 00 > . DV < 1 00 ) , DZ < 1 00 ) , CP <100) . T A < 650 ) 

EQU I VALENCE < KM , WS < 2 1 1 ) ) . < KM , MK ), <NP, WS < 20 1 ) ) . < SVM , WS < 2 1 0 ) ) 
1,<V2.Y3> 

INTEGER SVM ,P 

A. READ INPUT 
WRITE <6, 5) 

6 FORMAT <3 1 4 . 4F 12 . 5 ) 

8 FORMAT <3F 12.5) 

5 FORMAT <49H0XVZ POTENTIAL FLOW PROGRAM SECTION 6, 

? FORMAT < IX, QF 12.5) 

20 FORMAT < 1H1 , I4.31H STREAMLINES TO BE COMPUTED AT 

1 .F8.4.28H T FOR AN ONSET VELOCITY OF ,3F8.4) 

2 1 FORMAT < IX, 1 5HSTART I NG PC I NTS / , 3X , 2HPT , 5X , 1 HX , 1 1 X , 1 HV , 1 1X , 1 HZ ) 

22 FORMAT < IX, 14, 3F 12. 5) 

READ <03 ) <PF;OB< I ), ! = 1 . 15 ) 

WRITE < 6 . 90 > < PROB < I ), 1 = 1, 15) 

90 FORMAT < 1H0 , 15A4) 

B. READ THE PARAMETERS, T ARRAY AND SOURCE FROM TAPE 31 
READ <03 ) <WS< I ), 1=1,220) 

READ <03 ) <TX< I ), TV< I ), TZ< I ), XN< I ), VN< I >, ZN< I ),TA< I ), i=1, NP ) 

READ <03 ) SKIP 

IF < US <220) .EQ. 2. ) RE AD <03 ) SKIP 

READ <03 ) <S1< I >, 1=1 ,NP) 

READ <03 ) <S2< I ), l = 1^NP) 

READ <G3) <S3< I ), I = 1 ,NP ) 

C. READ THE B ARRRV 
WZ = Nr 

WZ=<WZ+1 1 .0/12.0 
NB = WZ 
IS = 2 
!F=24 1 

DO 12 IP = 1 , NB 
READ <04 ) F', <B< I 
I S= I S+240 
12 I F= I F+240 



AP< 1 )= .5 
AP<2) = .5 
AP<3) = 1 . 

AP <4 )=0 . 

AP<5)=0. 

GM< 1 ) = 1./6. 

GM<2) = 1./3. 

GM<3) = 1./3. 

83 READ <5,6 ) NOBP . NST , I END , DT , VX I , VY I , VZ I 
VSO=VX I *+2+VV I ++2+VZ I *+2 
NOB=NOBP 
DO 10 1=1 .NOB 

READ <5, 8 ) XP< I ), YF'< I ),ZP< I ) 



IF (E0F<5 ) .EQ. 0. ) GO TO 10 
NOBP= I - 1 

NR I TEC 6,9) NOBP, NOB 

0 FORMAT C 1HQ, 15, 28H STREANLINES SPECIFIED NOT , 13) 

GO TO 11 

10 CONTINUE 

11 CONTINUE 

NR I TEC 16 > NOBP, NST, I END , UX I , UV I , UZ I 
NR I TEC6, 20 > NOBF', NST, DT, UX I , UV I , UZ I 
NRITEC6/21 ) 

NFll TEC6.22) < I , XF'C I ), VPC I ), ZPC I ), 1 = 1, NOBP) 

C NOBP - NUMBER OF STREAMLINES TO BE TRACED. 

C NST - NUMBER OF STATIONS AT NHICH STREAMLINES SHOULD BE COMPUTED. 
DO 15 1=1, NOBP 
XTC I ) = XF'C I ) 

VTC I ) = VPC I ) 

ZTC I ) = ZPC I ) 

SKXCI >=0. 

SKVC I ) = 0. 

15 SKZCI ) = 0. 

I TC=0 
I RK=5 
98 K= 1 
P = 1 
J=1 

DO 100 1=1, NOBP 

C D. SET THE PARTIAL UELOCITV TO THE FREE STREAM UELOCITV 

UX 1 C I >=-1.0 
UV 1 C I >=0 . 

UZ1CI )=0. 

uxzci :>=o. 

UV2C 1 )=- 1 . 0 
UX3CI >=0.0 
UV3C I >=0.0 
UZ3C I )=- 1 . 0 
100 UZ2CI >=0. 

C E. START LOOP OUEFi THE OUADS. 

295 J=2 

C FI PICK UP QUAD. INFORMATION 

296 X 1=BCJ ) 

V 1=BC.J+ 1 > 

X2=BC J+2 > 

V2-BCJ+3) 

X3=BC J+4 ) 

X4=BCJ+5) 

V4=BCJ+6> 

XC=TXCK ) 

VC=TVCK> 

ZC=TZCK) 

A =TRCK > 

XX=BC J+7) 

VX=BCJ+8) 

2X=BCJ+9) 

XV-BCJ+ 10 ) 

VV=B< J+ 1 1 ) 

ZV=BCJ+ 12) 

C F2 COMPUTE LENGTH OF SIDES OF QUAD. 

D 12=SQ2F CX 1 , X2 , V 1 , V2 , 0 . , 0 . ) 

D23=SQ2F C X2 , X3 >2 , V3 , . 0 , . 0 > 

D34=SQ2FCX3,X4,V3,V4, .0, .0) 

D4 1=SQ2F CX4.X1.V4. V 1 , . 0 . . 0 > 

F3 COMPUTE SLOPE OF SIDES 



C 



I F < X2-X3 >305 , 300 , 305 

300 Cl 23=1. 

GO TO 310 

305 01123= < V2-V3 )/< X2-X3 ) 

Cl 23=0. 

310 I F < X3-X4 >3 15,311,315 

311 0134=1. 

GO TO 320 

315 CM34=<Y4-Y3 >/(X4-X3 ) 

0134=0. 

320 I F < X4-X 1 >325 ,321, 325 

321 0141=1. 

GO TO 330 

325 CH4 1 = < Y 1 -V4 ) / < X 1 -X4 > 

0141=0. 

330 !F<X1~X2 >335,331 .335 

331 Cl 12=1. 

GO TO 340 

335 OH 1 2 = < Y2- V 1 > / < X2-.X 1 > 

Cl 12=0. 

0 F4 COMPUTE QURDRRPQLE MOMENTS 

340 C I XX=B< J+ 17 ) 

0 1 XV=B< J+ 18 > 

0 1 YY=8< J+ 19 > 

0 F5 COMPUTE SIM AND COS OF SLOPE ANGLE FOR EACH. SI 

CV12=(Y2-Y1 >/D 12 
CV23= ( V3-V2 >/D23 
CV34=<V4-V3 >/D34 
CV41=<Y1-Y4)/D41 
CX 12-CX 1-.X2 )/D 12 
0X23= < X2-X3 > /D23 
CX34=(X3-X4 >/034 
CX41=<X4-X 0/041 

C F6 COMPUTE MAX D I AG INAL 

ST=S02F ( X 1 , X3 , VI, Y3,0. ,0. ) 

8T2=SQ2F < X2 , X4 , V2 . V4 , . 0 , . 0 > 

I FCST-ST2 >341,342,342 

341 ST=ST2 

C: G. START LOOP OYER THE OFF BODY POINTS 

342 DO 530 JQ=1,NQBF 
I S= 1 

XCQ=XP<JQ> 

YCQ=VP<JQ> 

ZCQ=ZF<JQ) 

J 1= 1 

345 RPQ=SQ2F ( XC . XCQ . YC . VCQ , ZC , ZCQ > 

C H. DETERM IN METHOD 

! F(RP0-ST*4 >350, 350 . 450 

350 X= ( XCQ-XC )*XX+ C VCQ- YC >*VX+ C ZCQ-ZC )*ZX 
V= C XCQ-XC >*XY+ «: VCQ- YC >*YY + < ZCQ-ZC )*ZY 
Z=< XCQ-XC >*XM(K >+(YCO-YC )*VN<K )+( ZCQ-ZC )*ZN<K > 

1 F<RPQ-ST*2 . 5 >355 , 355, 400 

C I. COMPUTE I MOUSED VELOCITY BY EXACT METHOD 

355 R 1 =SQ2F < X , X 1 , Y , V 1 , Z , 0 . > 

R2=SQ2F O' , X2 , V , V2 , Z , 0 . > 

R3=SQ2F (X, X3, V, Y3, Z, 0 . > 

R4=SQ2F <X , X4 , Y , Y4 . Z . 0 . > 

CLA 1=AL00< <R 1+R2-D 12 >/<R 1+R2+D 12 > > 

CLA2=AL0G< (R2+R3-D23 >/<R2+R3+D23 > > 

CLA3=RL0G< CR3+R4-D34 >/(R3+P.4+D34 > > 

CLR4=RL0G< <R4+R 1-D4 1 >/<R4+R1+D4 1 > > 

TUX=CV 12+CLA 1+CV23*CLR2+CY34*CLR3+CV4 1+CLA4 
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TU'r-CX 1 2*CLfi 1 +CX23*CLR2+CX34*CLR3+CX4 1*CLR4 
TUZ=U . 

I F < RBS < Z >- . 00 1 *ST >375,361,361 
361 ZSCl=Z*+2 

E 1=ZSQ+<X-X 1 >**2 
E2=ZSQ+(X-X2 )**2 
E3=ZSQ+ <X-X3 >**2 
E4=ZS0+(X-X4 )**2 

hi=<v-v n*a-xn 

H2=<V-V2 >*<X-X2> 

H3=<V-V3)*«-X3> 

H4=<V-V4 )*<X-X4 > 

! F < C I 12)363,363,364 

363 US 1 = < CH 1 2+E 1-H 1 )/( Z*R 1 > 

MS2= < CM 1 2*E2-H2 >/<Z*R2 ) 

AT 1=ATRN(MS1 > 
fiT2=RTRN<US2> 

TUZ=RT 1 -RT2 

364 I FCC 1 23 >366,366, 367 

366 RT3=RTRN < ( CM23+E2-H2 > / < Z*R2 ) > 
ftT4=RTF:MC <Ct123*E3-H3 >/<Z+R3 > > 

TUZ=TUZ+fiT3-ftT4 

367 1 FCC i 34 >368 , 368 , 368 

363 AT5=RTRriC <C1134*E3-H3 >/<Z*R3 ) ) 

P.T6=fiTfiN«Ci134*E4-H4 >/<2*R4 > ) 

TUZ=TUZ+RT5-AT6 

368 !F(C 141)370,370,375 

370 RT7=RTRN < < 0114 1 +E4-H4 > / < Z+R4 > > 

RT8=RTRM( <CM4 1*E 1-H 1 >/<Z*R 1 > > 

TUZ=TUZ+RT7-RT8 
375 GO TO 450 

C J. COMPUTE I HOUSED UELOCITV BV QURDRRPOLE METHOD 

400 RPQ3=RPQ**3 

RPQ7= < RPQ3**2 >*RPQ 
US1= X/RPQ3 
XSQ=X**2 
VSQ=V**2 
2SQ=Z* :t: 2 

PS=VSQ+ZSQ-4 . *XSQ 
QS=XSQ+ZS0-4 . *VSQ 
WS2=X* < 9 . *PS+3Q . *XSQ > /RPQ7 
US3=3 . *V*PS/RPQ7 
WS4=3 . *X*QS/RPQ7 

TUX=fl*US 1 -C I XV+US3-C I XX*US2-C I VV*US4 
US 1=V/RP03 

US2=V*<9 . *03+30 . *YSQ J/RPQ7 
TW=fl*US 1-C I XX+WS3-C I XV+U34-C I VV+MS2 
TUZ=Z* ( R /RPQ3-3 . *<C I XX*PS-5 . *C I XV*X*V+C I YY*QS >/RP07 > 
450 UXC I S >=TUX*XX+TUV+XV+TUZ*XN < K ) 

UV< I S >=TUX*VX+TUV*VV+TUZ*VH < K > 

UZ< I S >=TUX*ZX+TUV+ZV+TUZ+Z1KK ) 

GO TO 470 

C K. COMPUTE I MOUSED UELOCITV BV MOMOPOLE METHOD 

460 RRPQ3=fi/<RPQ**3> 

UX < I S >= < XCO-XC >+RRF'Q3 

W(IS>= < VCQ-VC >*RRPQ3 
UZ< IS>= CZCQ-ZC )+RRPQ3 

C L. REFLECT OFF BOOV POIMT IH PLANE OF SVMETRV 

470 GO TO < 480 , 485 , 490 , 495 , 500 , 505 , 5 1 0 , 5 1 5 > , I S 
480 UKJ1>=UXa> 

u 2 (Ji>=uxa> 

U3(U1)=UX(1> 



UKJ1+1 )=UY< 1 ) 

U2<J1+1 >=UY< 1 ) 

U3<J 1+ 1 )=W< 1 ) 

UKJ1+2)=UZ< 1 ) 

U2< J 1+2 >=UZ< 1 ) 

U3(J1+2>=UZ< 1 ) 

IF(SVM) 525,525,481 

481 18=2 

C XZ SVMETRV 

VCQ—VCQ 
GO TO 345 

4S5 I FCSYtl- 1)517,517,486 
C XV SVMETRV 

486 18=3 
ZCQ— ZCQ 
GO TO 345 

490 18=4 

VCQ=-VCQ 
GO TO 345 

495 I FCSVM-2 >516,516, 496 

C V2 SVMETRV 

496 18=5 
XCQ=-XCQ 
GO TO 345 

500 18=5 

VCQ=-VCQ 
GO TO 345 

505 18=7 

ZCQ=-ZCQ 
GO TO 345 

510 18=8 
VCQ=-VCO 
GO TO 345 

C M. ROD CONTRIBUTIONS OF RLL REFLECTIONS 

5 15 U 1 < J 1 )=U 1 ( J 1 )+UX ( 8 >+UX(7 >+UX(6 >+UXC5 ) 

U2< J 1 >=U2< J 1 >-UX<8 HUX<7 >+UX<6 )-UX<5 ) 

U3 ( J 1 >=U3 < J 1 >-UX < 8 )-UX < 7 >+UX < 6 >+UX ( 5 ) 

UK J 1+1 >=U K J 1 + 1 >-UV ( 8 )+U V (7 >+UV < 6 ;-UV < 5 ) 
U2<J 1+1 >=U2< J 1+ 1 >+UV<8 >+UV<7 >+W<6 >+UV<5) 
U3C J 1+ 1 )=U3(J 1+ 1 J+UVC8 ?-UV < 7 HUY ( 6 )-UY (5? 

U 1 ( J 1+2 )=U 1 ( J 1+2 )-UZ(S )-U2<7 )+U2<6 )+U2<5 ) 
U2C J 1+2 )=U2(.J 1+2 >+UZ<8 )-U2C7 >+UZ<6 )-U2<5 > 
U3 < J 1 +2 )=U3 ( J 1 +2 )+UZ ( 8 >+U2 < 7 )+U2 < 6 )+U2 (5 ) 

516 U K J 1 )=U 1 < J 1 >+UX<4 >+UX(3 ) 

U2<U 1 )=U2< J 1 )+UX(4 )-UX ( 3 ) 

U3CJ 1 >=U3CJ 1 >-UX(4 >-UX(3 > 

U K J 1 + 1 >=U K J 1 + 1 >+UY < 4 >-UV < 3 > 

U2< J 1+ 1 )=U2< J 1+ 1 ?+UVC4 >+UY (3 ) 

U3CJ1+1 )=U3( J 1+ 1 >-UV<4 >+UV<3 > 

U K J 1 +2 >=U K J 1 +2 >-UZ ( 4 >-U2 < 3 ) 

U2 C J 1+2 >=U2 < J 1 +2 >-UZ < 4 >+UZ < 3 > 

U3 < J 1 +2 )=U3 < J 1+2 )+U2<4 >+U2 < 3 ) 

517 U K J 1 >=U K J 1 HUX < 2 ) 

U2(J1>=U2(J1HUX<2> 

U3 < J 1 >=U3 < J 1 >+UX < 2 ) 

UKJ1+1)=UKJ1+1)-UVC2) 

U2CJ1+ 1 >=U2< J 1+ 1 HUY <2 ) 

U3CJ1+1 >=U3(J1+1 >-UY<2> 

U K J 1+2 )=U 1 < J 1+2 HUZC2 ) 

U2 < J 1 +2 >=U2 < J 1 +2 )-U2 < 2 ) 

U3 < J 1 +2 )=U3 < J 1 +2 )+UZ < 2 > 

525 UX 1 ( JQ HUX 1 < JQ )+U 1 < 1)*S 1 < P > 
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W 1 C JQ >=UV 1 < JQ )+U 1 <2 )+S 1 CP ) 

UZ 1 ( JO >=UZ 1 < JQ >+U1(3>+S1(P> 

UX2 < JQ >=UX2 ( JQ >+U2 C 1 )+S2 < P ) 

UV2 C JQ )=UV2 < JQ >+U2 C 2 >+S2 < P ) 

UZ2 < JO >=UZ2 < JQ >+U2 < 3 )+S2 < p > 

UX3< JQ >=UX3( JQ >+U3< 1 )*S3<P ) 

UV3 ( JQ >=UV3 < JQ >+03 ( 2 >+83 < P > 

530 UZ3 C JQ >=UZ3 C JQ >+U3 C 3 >+S3 C P > 

C N. END OF LOOP OUER OFF BuDV POINTS 

585 P=P+1 
K=K+ 1 
J=J+20 

I F < K-NP >295 , 296 , 599 
C P. END OF LOOP' OUER QUADS 

599 H=fiP< IRK >+DT 

DO 730 I = 1 , NOBP 
63 FORMAT C2X , I3,3F12.5,9X,4F12.5> 

DX ( I >=- (UX I +UX 1 C I >+UV I +UX2 ( I >+UZ I +UX3( I > > 

DV C I >=- < UX I +UV 1 C I >+UV I +UV2 C I >+UZ I +UV3 < I >> 

730 DZ< I >=-<UX I +UZ 1 < I )+UV I +UZ2 ( I >+UZ I +UZ3C I > > 

I F < I RK . EO . 5 > GO TO 900 

I F ( I RK . EQ . 4 > GU T 0 SOU 

DO 750 1=1, MOBP 

XP< I >=XT < I )+DXC I >+H 

VPC I >=VT< I )+DYC I )+H 

ZP(I >=ZT < I >+DZ< I >*H 

SKX< I >-8KX< I >+GM< IRK>+DXC I > 

SKV< I )=SKV( I >+GM< IRK)*DV( I > 

750 SKZC I >=SKZ< I >+GH< IRK>*DZ< I > 

IRK = IRK + 1 
GO TO 98 
300 H=DT 

DO 330 1=1, NuBP 

DX < I >=- <UX I +UX 1 ( I >+UV I +UX2 < I >+UZ ! +UX3 ( I > > 

DV< I >=-<UX i +UV 1 < I >+UY I +UV2 ( I >+UZ I +UV3( I > > 

DZ ( I >=- ( UX I +UZ 1 < I >+UV I +UZ2 ( I >+UZ I +UZ3 < I > > 

XP< I >=XT< I )+H+DX ( I >/6 . +SKXC I >+H 
XT ( I )=XP< I > 

VP X I >=VT ( I >+H+DV( I >/6 . +SKV< I >+H 
VT< I >=VP< I > 

2P< I >=ZT < I >+H+DZ( I >/6.+SKZ< I )+H 
ZT< I >=ZP< I > 

SKX( I >=0. 

SKY< I ) = 0. 

830 SKZ < I > = 0. 

IRK = 5 
GO TO 98 
900 IRK = 1 

DO 905 1=1, NOBP 

DSQ=OXC I )++2+DV< I >*+2+DZ( I >++2 
CPC I >= 1 . -030 /USQ 
905 CONTINUE 

URITEC6,61 > 1 TC 

61 FORMAT <6HQ STEP. 1 4/ > 

URITEC6,62> 

62 FORMAT C3X.4HL I NE.5X, 1HX, 1 IX , 1HV, 1 IX, 1HZ, 20X, 2HUX, 10X,2HUY. 

1 10X, 2HUZ , 10X, 2HCP) 

UR I TEC 6,63) < I ,XP< I >, VP< I >, ZP< I >, DX< I > . DVC I >, DZC I > . CPC 1 >, 1 = 1, HOBP > 
UR I TEC 16 > < XP< I >, YP< I > , ZF'C I >, 1 = 1, NOBP > 

I F < ITC.EQ.NST) GO TO 910 
I TC= I TC+ 1 
GO TO 599 
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910 I F c I END . EQ . 0 > GO TO 83 
REWIND 03 
REWIND 04 
ENDFILE 16 
REWIND 16 
STOP 6 
END 

FUNCT I ON SQ2F (X 1 , X2, V 1 , V2 . 2 1 , Z2 ) 
X=X1-X2 
V=V1-V2 
2 = 21-22 

RS=2**2+V**2+X**2 

F;=RBS(X )+RBS< V >+RBS<Z >+ 1 . 0E-20 

R= . 3422* < R+ < RS+FlS )/R ) 

R= R+RS/R 

SQ2F= . 25*R+RS/R 

RETURN 

END 
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APPEND I X VII ~ XVZPF SECT I OH PF7 



PROGRAM PFP7<TRPE7, I HPUT= 128, OUTPUT* 128, TfiPE5= I NPUT, TfiPE 17, 
1TRPE6=0UTPUT, TRPE03, TRPE04, TRPE3=TflPE03, TfiPE4=TfiPE04 > 

XVZ POTENTIAL FLOW PROGRAM UERSION 4 AND HERS I ON 5 SECTION 7 
COMPUTES VELOCITIES AMD PRESSURE COEFFICIENTS FOR 
ON BODY STREAMLINES 



COMMON X < 658 > , Y ( 658 > , Z < 658 ) , XH < 650 ) , YN <650 > 

1 , ZN<650 ), UX 1 <650 ), VV 1 <650 ), UZ 1 <650 ), VX2C650 ), VY2<650 > 

2, UZ2< 650 ) , UX3 < 650 ) , UY3 < 650 ) , UZ3 < 650 > , XC 1 < 658 > , VC 1 < 658 > 

3 , XC2 < 658 ) , YC2 < 658 ) . XC3 (658), XC4 < 658 ) , YC4 < 658 ) , 

4X3 < 658 ) , V3 < 658 > , Z3 < 658 ) , X4 < 658 ) , Y4 < 658 ),Z4<658) 

D I MENS I ON XL < 1 50 ) . YL < 1 50 ) , ZL < 1 50 > , UX < 150 > , UY < 1 50 > , UZ < 150), 

1CPC 150 ), GK 1 < 150 ) , GK2< 150 ), H2< 150 ), STML< 150 ), UASS< 150), NOUAD< 150 ) 
5 , DMX < 650 ) . F'FiOB < 1 5 > . YC3 < 658 ) , SF < 5 > . XCR < 5 ) , VCR < 5 ) 

7HSP < 50 ) , US < 220 ) , XST < 50 ), YST < 50 ) , ZST < 50 ) ' 

EQU I VALENCE < US <201), NP ) , < YC3 , VC2 ) 

RERQ<G3XPRGB< I ), 1 = 1, 15) 

READ <03 XUS< I ), 1=1, 220 ) 

READ <03 )<X< I ), V< I ),Z< I ), XN< I ), VN< I ) . ZN< I ) , 

1SKIP. 1=1, NP) 

READ<03)SKIP 
IUER =US<220) 

WRITE<6,5) IUER 
5 FORMAT < 4 6H0XVZ POTENTIAL 
IF < I VEFi . EQ . 5 ) RE AD < 03) 

RERD<03)SKIP 
READ < 03 )SK I P 
RERD<03)SK IP 



FLOW PROGRAM SECTION 
SKIP 



VERSION 



READ <03 XVX 1< I ),VY1< I ), VZX I ), 1 = 1, NP ) 
READ <03 X VX2< I ) , UY2< I ) , VZ2< i ) , I = 1 , NP ) 
READ <03 XVX3C I )‘ VY3< I 'i , VZ3< i >’ I = 1, NP ) 
REWIND 03 
NB=<HP+ 11)/ 12 
DO 80 1=1, NB 
I FN= I + 12 
1 8= ! FN- 1 1 



READ <04 ) Q, < XC 1 < J ) , VC 1 < J ) , XC2 < J ) , 

1 YC2 < J ) , XC3 < J ) , XC4 < J ) , YC4 < J > , X3 < J ) , Y3 < J ) , 

2Z3 < J ) , X4 < J > , V4 < J > , Z4 < J > , <8KIP,K=1,7), J=IS, I FN ) 
NQ=Q 

IF<NQ.NE. IS) GO TO 450 
80 CONTINUE 
REWIND 04 
DO 90 1=1, NP 

D 1=<XC 1 < I )**2+VC 1 < I )*" t: 2 X 1.01 
D2=<XC2< I )**2+VC2< I )**2 )+ 1 . 0 1 
D3=<XC3< I )**2+YC3( I )**2 )* 1.01 
D4=<XC4< I >**2+VC4< i >**2 >* 1.01 
90 DMX< I )=RMRX 1 <D 1 , 02, 03 , 04 ) 

1 1 FORMAT <3F 12 . 4, 3 1 4, F 12 . 4 ) 

12 FORMAT <3F 12 . 4, 14) 

M 1 0=75 



100 READ<5, 1 1 ) UX I . VV I . VZI , NL IN,MRXJ. IWRITE.RMACH 
IF <EOF<5).NE.O. ) NL I N=Q . 

HXJ-NAXJ 

IF < MAXJ .LE. 0 .OR. MAXJ .GT. NP/2) MAXJ = NP/2 
M I NJ=M I O-MAXJ 
MAXJ=H I D+MAXJ 

IF <MAXJ . GT . M 1 0*2 ) MAXJ=M 1 0*2 



IF < MINJ LT. 1) MIHJ = 1 
WRITE <7 > NUN 
WRITE-,' 17) NLIN,UXI,UVI,UZI 
IFCNLIN.LE.O) GO TO 550 



WRITE<5,30> <PR0B( I ), 1=1, 15) 

WR I TE ( 6 ‘ 34 > UX I , W I . UZ I .HUH, NX J , I WR I TE . RfiftCH 
34 FORMAT <34H00N BODV STREAMLINES - INPUT DATA /6H UX I =,F10.5/ 
16H UVI =, F 10 . 5/OH UZ I =,F10.5/5H NLIN=, I 10/ 

2 6H JMAX=, I 10./.SH IWRITE=, ! 10, /,QH MACH N0=,F10.5> 
URITE<6,38> 



FORMAT (27H0STFiEfif1L I ME STARTING P0INTS/5H LINE, 1 1X, 1HX, 12X, 1HV, 



1 12X, 1H2, 10X,3HNSP> 



L I N=HL I N 



00 45 1=1, LIN 

READ <5, 12) XST< I ) ,VST< I >,ZST< I >.NSP< I > 
IF CE0F(5).E0.0. ) GO TO 45 
NL I N= I - 1 



WRITE <6, 42) NUN, LIN 

42 FORMAT < 1 HO , 1 5 , 28H STREANL I NES SPEC I F I ED NOT ,13) 
IF(NLIN.LE.O) GO TO 550 



GO TO 48 

45 WRITE <6, 45) I .XSTC I ), VST( I J.ZSTC i ),NSP< I > 
45 FORMAT < IX , 1 3, 2X , 3F 13.5, 19) 

48 CONTINUE 

USQ=UX I ++2+UV I **2+UZ I **2 
IF (fiMACH .EQ. 0. ) GO TO 1130 
C *** COMPUTE CRITICAL MACH NO. 

USD = 0. 



DO 1100 1 = 1, MF* 

US = ( UX I +UX 1 < I )+UV I *UX2 < I )+UZ I *UX3 < I ) )**2 + 

1 < UX I *UV 1 < I )+W I *UV2 < I )+UZ I *UV3 < I > )**2 + 

2 ( UX I *UZ 1 < ! )+UV I :f UZ2< I >+UZ ! *UZ3 < I ) )**2 
IF (US .GT. USD) USD = US 

1100 CONTINUE 

U = SORT (USD /USQ) 

CMNA = 1 . /U 
DO 1110 1=1,3 

CMNB = <<<CMNA**2+5. )/6. >++1.75)/U 
CMNC = < < < CMNB+*2+5 . ) /5 . >** 1 . 75 ) /U 
1110 CMNA = <CMNR*CMNC-CMN8**2 >/<CMNfl+CMNC-2 . *CMN£ ) 
WRITE-: 5, 1 120) CMNA 

1120 FORMAT <2 1H CRITICAL MACH NO. =,F5.3) 

1130 CONTINUE 

C START LOOP OUER STREAMLINES 
DO 400 LL=1.MLIN 
D I RT= 1 . 

101 Jl=1 
RF= 1 . 

UX(MID)=0. 

IJV(MID)=0. 

UZ (MID )=0 . 

CP < M I D )=0 . 

H2(M I D )= 1 . 

GK 1 CM ! D >=0 . 

GK2(MID)=0. 

STML<MID)=0. 

102 NQ=NSP<LL> 

LNQ=NQ 

xl(mid)=xst-;ll) 

VL<MID)=VST<LL) 

ZL(MI D)=ZST (LL ) 

J=MID 
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•JL=J 

C SEPARATE CALCULATION OF SECOND 

C POINT FROM MAIN LOOP 

XLT=(XL< J )-X<NQ ) )*X3(NQ )+(VL (J >-Y(NQ ) )*V3<NQ ) 

1 + < ZL ( J )-Z ( NQ ) )*Z3 < NO ) 

YLT= C XL < J )-X < NQ ) )*X4 < NQ >+ < YL < J )-Y < NQ ) )*Y4 < HQ ) 

1 + < 2L < J )-Z < NQ ) )*Z4 < NQ > 

XL< J )=XLT*X3(NQ H7LT*X4(NQ HX(NQ > 

YL< J )=XLT*Y3<NQ )+VLT*Y4 (NQ HY(NQ ) 

ZL< J )=XLT*Z3(NQ HYLT*Z4(NQ >+Z(NQ > 

105 IOT=NOD(NQ, 4 ) + 1 

GO TO <630,600,610,620) IQT 
600 NR-NQ+I 
NU=NQ+2 
GO TO 10? 

610 NR=N0+2 
NU=NQ- 1 
GO TO 107 
620 NR=N0-2 
NU=NQ+ 1 
GO TO 107 
630 NF;=fTQ- 1 
NU=NQ-2 

1 07 UXQ=- < UX I *UX 1 < NQ HUY I *0X2 < NQ HUZ I *0X3 < NQ > > 

UYQ=-<OX I *UY 1 (NQ HUY I *UY2(NQ >+UZ I *073<NQ ) > 

UZO=-<OX I *UZ 1 (NQ )+UY I *UZ2(NQ >+0Z I *023 (NQ > > 

UXFl=-(OX I *0X 1 (NR )+0V I *0X2 (NR HUZ 1 *0X3 CNF; > ) 

UYR=- ( OX I *0V 1 ( NR HUY I *0V2 ( NP> )+0Z I *UV3 < NR ) ) 

UZF;=-(OX I *0Z 1 (NR HUY I *0Z2(NFi HUZ I *0Z3(NR ) ) 

UXU=-(OX I *UX KNU HUY I *0X2 (NU HUZ I *0X3 (NU ) ) 

UVU=-(OX I *0Y 1 (NU HUY I *0V2(NU HUZ I *073 (NU ) ) 

UZ!J=- ( OX I *02 1 ( NU HUY ! *022 < NU >+0Z I *UZ3 ( NU ) ) 

C TRANSFORM UELOCITIES TO QUAD SYSTEM 
UQ=UXQ*X3(NQ HUYQ*'Y3(NQ >+U2Q*Z3(NQ > 

UQ=UXC*X4(NQ HUVQ*Y4(NQ )+UZQ*Z4(NQ ) 

CSR= 1 . /<XN(NQ )*XN(NF; )+VN(NQ)*VN(NF; HZN(NQ )*ZN(NR ) ) 
UT=UXR*X3(NR HUYR*Y3(NR )+UZR*Z3(NFl ) 

UT= ( UXR*X4 ( NR )+UVF;*'V4 < NR )+UZR*Z4 (NR) )*CSR 
XXFi= < X3 ( NR )*X3 < NQ )+ V3 < NR )*V3 ( NO HZ3 < NR >*Z3< NQ ) > 

XVR= ( X4 ( NR )*X3 ( NQ )+ V4 ( NR )*V3 ( NQ HZ4 ( NR )*Z3 < NQ ) ) 

UR=UT*XXR+UT*XYR 

YXR= < X3 < NR )*X4 < NQ HY3 ( NR )+Y4 < NQ HZ3 < HR >*24 < NQ ) > 

VVF;= ( X4 < NR )*X4 ( NQ )+ V4 ( NR )*V4 ( NQ HZ4 < NR )*Z4 ( NQ ) ) 

OR=UT*VXR+OT*VVR 

UU=UXU*X3(NQ HUVU*V3(NQ )+UZU*Z3(NQ ) 

CSU= <XN(NQ )*XN(NU )+VN(NQ )*YN(NU HZN(NQ )*ZN(NU ) > 
0U= ( UXU*X4 ( NQ HU YU*Y4 < NO >+UZU*Z4 ( NQ ) ) / CSU 
C FIND RELAT10E COORD INRTES OF NEIGHBORING OUADS 
XD=X(NF;)-X<NQ) 

VD=V(NR HY(NQ) 

ZD=Z(NFl)-Z(NQ) 

XT =XD*X3 ( HR )+ VD+V3 ( NR )+ZD*Z3 ( NR ) 

YTT =XD*X4 ( NF; )+YD*V4 ( MR )+ZD*Z4 < NR > 

ZT =XD*XTi < NR )+ VD*YN C NR HZD*ZN < HR ) 

YT=(-4*SQRT(VTT**2+ZT**2 )+VTT*CSR+YTT )*CSR* . 16666667 

XR=XT*XXP.+YT*XVR 

VR=XT*YXR+VT*VVR 

XD=X(NIJ)-X(NQ) 

VD=Y(NU )-V(NO) 

ZD=Z(NU )-Z(NQ > 

XU=XD*X3(NQ )+YD*V3(HQ )+ZD*Z3(NQ ) 

VT=XD*X4 ( NQ )+YD*Y4 ( NQ )+ZD*Z4 < NQ ) 
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2T=XD*XN < NQ )+YD*YN < NQ )+ZD*ZN< NQ > 

VU=<4 . *SQRT<YT**2+ZT**2 )+VT/CSU+VT )* . 1666656? 

C FIND COEFFICIENTS OF VELOCITY FUNCTIONS 
DEN=1 . /<XR*VU-XU*VR) 

U 1= < < UR-UO )*YU- ( UU-UQ >*YR > ♦DEN 
U2=- < ( UR-UQ )*XU- < UU-UQ >+XR >*DEN 
U 1=( <UR-UQ >*VU-< UU-UQ >*YR ) *DEN 
U2=-( ( UR-UQ )*XU- < UU-UQ )*XR )*OEN 
C FIND UELQC I TV RT STREAMLINE POINT 
USL»UQ+U1*XLT+U2*YLT 
USL=UQ+U 1+XLT+U2+VLT 
UXP=USL*X3<NQ )+USL*X4(N0 ) 

UYP=USL*Y3<NQ >+USL+V4(NQ > 

UZP=USL*Z3<NQ >+VSL*Z4(NQ > 

C FIND GEODESIC CURVATURES GK 1 , GK2 
VS0D=USL**2+VSL**2 
DEN=USQD*SQRT<VSQD ) 

GK 1P=CUSL*<VSL*U2-USL*U2 HJSL*<USL*U 1-USL*V 1 > )/DEN 
GK2F- < USL* < USL*U 1 +USL+U2 >-USL* < USL*U 1 +VSL*V2 ) > /DEN 
C FIND LOCAL STREAM FUNCTION 

CXV-CU 1 *00**2- v2*UQ ++ 2 > /USQD 
CVY*U2-UQ*VQ* < U1+U2) /USQD 
CXX=IJ2-CVV-U 1 

C0=XLT*U0-VLT*U0-CXV+XLT*VLT-CVV*VLT+*2-CXX*XLT**2 
C FIND STREAM FUNCTION AT CORNER POINTS 
XCRC1)-XC1<NQ> 

XCFK2 )=XC2<NQ) 

XCR<3 )=XC3<NQ ) 

XCR < 4 >=XC4 < NQ > 

XCR<5 )=XCR< 1 > 

VCR( 1 )=VC 1<NO > 

VCR C 2 >=VC2 < NQ ) 

YCF;C3)=YC3CNQ) 

VCR (4 )=VC4(HQ> 

VCR<5>=VCR< 1 > 

DO 110 N=1,4 

1 1 0 SF < N >=CO-VQ*XCR < N )+UO*VCR < N >+CXV*XCR < N >*VCR ( N >+CYY*YCR < N >♦•* 
1 +CXX*XCR<N>**2 
SF<5>=SF< 1 > 

TEST=0 

DO 120 N=1,4 

IF < SF(N)*SF<N+ 1 ) .GE. 0. > GO TO 120 
Xf1=(XCFi<N HXCFKN+ 1 ) >* . 5 
VM=CVCRCN )+VCR<N+ 1 > )* . 5 
C FIND INTERSECTION WITH SIDE OF QUAD. 

SFM=C0-VQ*XM+UQ*YM+CXY*Xf1*YM+CYY*YM**2+CXX*XM**2 
RC=2 . * < SF< N )-2 . ♦SFM+SF < N+ 1 ) > 

BC=3 . *SF ( N )-4 . ♦SFM+SF <N+ 1 > 

IF (AC .EQ. 0) GO TO 112 
SR=SQRT< BC**2-4 . *AC*SF <N> ) 

TF- < BC+SR > / < 2 . ♦AC ) 

IF <TP .LE. 1. .AND. TF .GE. 0. ) GO TO 115 
TF- < BC-SR > / < 2 . ♦AC ) 

GO TO 115 

113 IF <BC .EQ. 0 ) GO TO 120 
TP=SF(N)/BC 

1 15 XNF'=< 1 . -TP >*XCR<N )+TP*XCRCN+ 1 > 

VHP=C 1 . -TP >*YCR(N >+TP*YCR<N+ 1 > 

TESTP= ( < XNP-XLT >*UQ+ < VHP-VLT >*UQ >♦0 1 RT 
IF (TESTP .LE. TEST) GO TO 120 
TEST=TESTP 
XNT =XNP 
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YNT=YNP 
120 CONTINUE 

I F< TEST .EQ. 0) GO TO 280 
C RUERRGE LRST UELOCITV RND CURURTUP.E 
UXCJ>=<UX<J)+UXP)*fiF 
UV(J>=<UY<J>+UYP>*flF 
UZ < J >= < U2 < J )+U2P >*RF 
GK K J >= < GK 1 < J >+GK 1 F >*RF 
GK2 < J )= < G'K2 < J >+GK2P >*RF 

H2< J >=H2< JL )*<2 . -GK 1 < JL >*($TML<J >-STML< JL > > )/<2 . +GK 1 < J )*< 
1 STML C J >-STNL < JL > ) > 

CP < J >= 1 . - CUX < J )+*2+llV ( J >**2+UZ < J >**2 > /USQ 
URBSC J >=SQF;T < 1 . -CPC J) ) 

C CONFUTE UELOCITV RT NEXT POINT 

NQURDC J >=NQ 
JL=J 
J=J+JI 

USL-UQ+XNT*U 1 +VNT+U2 
USL=UQ+XNT*U 1 +VNT+U2 
UX<J)=USL*X3CNQ )+USL*X4CNQ > 

UV < J MJSL*V3 < NO >+USL*Y4 ( NQ ) 

U2 < J )=USL*Z3 < NQ >+U3L*Z4 < NQ > 

C COMPUTE GEODESIC CURURTURES 

USQD=USL**2+USL**2 
DEN=USQD*SQRT CUSOD > 

GK 1 (J >= < USL* < USL*U2-USL*U2 >-USL*<USL*lJ HJSL*U 1 > ) /DEN 
GK2CJ )=<USL*<USL*L! 1+USL*U2 >-U$L*<USL*U 1*'JSL*U2 > >/DEN 
eGRD=SQF;T<<XNT-XLT)*+2+<VNT-VLT>**2 )*DIRT 
STML C J >=ST NL < JL )+CGF;D 
C COMPUTE H2 

H2 < J )=H2 < JL >* ( 2 . -CORD*GK 1 C JL » / < 2 . +CORD+GK 1 < J » 

CP<J)=1.-US0D/USQ 

URBSC J )=SQRT < 1 . -CP C J ) > 

RF= .5 
LNQ=NQ 

XL C J >= XNT *X.3 < N'Q )+VNT*X4CNQ >+XCNQ ) 

VL < J )=XNT*Y3 < NO >+VNT*Y4 C NO >+Y C NO ) 

ZL < J )=XNT*Z3 < NQ >+YNT+Z4 < NQ )+Z C NQ ) 

IF ( J .LE. Ml HJ .OR. J .GE. MfiXJ > GO TO 280 
C FIND NEXT QUAD. 

1 = 1 

250 HQ= I 

!F< I .ECi.LNQ) GO TO 280 
TEST=«L<J>-X< I > )**2+<VL< J )-V< I >>**2+ 

1CZLCJ)-ZC I )>*+2-Dt1X< I ) 

IFCTEST.GT.O. > GO TO 280 
DS1=CXC1C I >-XC2< I )>*+2+<VC1< I >-VC2C I »*+2 
D-S2-CXC2C I >-XC3< I »**2+CVC2< I )-YC3C I »**2 
DS3=CXC3C I )-XC4 < I >>**2+<YC3< I >-YC4C I >>**2 
DS4=(XC4< I >-XC 1 < I >>**2+<VC4< I >-VC1< I >)+*2 
XLT=CXL(J)-X( I >}+X3< I >+<VL<J)-V< I >>*V3< I >+ 

1CZLCJ )-ZC I >>*Z3CI > 

YLT=(XL( J )-XC I )>*X4C I )+CYLCJ)-VC I )>*V4C I )+ 

1 CZLCJ )-Z< I ) >*Z4C I > 

ZLT=<XL<J)-X< I )>*XN< I >+<YL(J >-V< I >)*VN< I )+ 

1 < ZL < J )-Z < I ) )*ZN C I ) 

ZSQ=2LT**2 
TE8T=ZS0- . 1+DMXC I ) 

IFCTEST.GT.O. > GO TO 280 

FiC 1 =SQRT C ZSQ+ < XLT-XC K I ) >**2+<YLT-VC K I > > 

1 ** 2 ) 

RC2=SQF;T < ZSQ+ < XLT-XC2 < I )>**2+<VLT-YC2< I ))+*2) 
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RC3=SQRT ( ZSQ+ ( XLT-XC3 C I )>**2+<VLT-VC3C 1 >5**2) 

RC4=SQRT CZSQ+CXLT-XC4C I ) )**2+<VLT-VC4< I ))**2 ) 

TEST= C CRC 1 +RC2 )**2 >-DS 1 *1.21 
IFCTEST.LT.O. ) 60 TO 105 
TEST = < < RC2+RC3 )**2 >-DS2 *1.21 
IFCTEST.LT.O. > GO TO 105 
TEST= C <RC3+RC4 >**2 >-DS3 *1.21 
IF<TEST.LT.G. ) GO TO 105 
TEST= < (RC4+RC 1 )**2 )-DS4 *1.21 
IFCTEST.LT.O. > GO TO 105 
280 1=1+1 

IFCI.LE.NP) GO TO 250 
282 IF (DIRT .LT. 0.) GO TO 285 
0 i RT=- 1 . 

J I =- 1 
JMRX=J 
GO TO 102 
285 JM I M=J 

SS=STNLCUNIN> 

DO 290 J=JMIN, JMRX 
290 STML ( J >=STML C J >-SS 
JHN=Jt1IN+1 
JMX=JMAX-2 
RF= 1 . 

L=.JMN 

URITEC6.30XPR0BCI >, 1 = 1, 15) 

30 FORMRT C 1 H 1 , 15R4) 

UR I TE (6 ,20 >0X1,071,021 

20 FOPiflRTC 18H0 ONSET FLOW, 0XI = ,F6.3,2X.4H0VI=,F&.3,2X 1 4H0ZI = ,F5.3) 
UR I TE < 6 . 50 > LL , NSP (LL), XST (LL ), VST <LL>, ZST <LL) 

50 FORfIRT < 1 1H0 LINE NO. , I2,31H PASSING THROUGH QUADRILATERAL ,13, 

1 28H WITH STARTING POINT, X=,F12.5,2X.2HV=,F12.5,2X, 

2 2HZ= ,F 12.5 //) 

IF CJMIH.LE.MINJ .OR. UMAX . GE . MAXJ > UR I TEC 5 , 55 > 

55 FORMAT (35H PROBABLE ERROR - LINE IS UERV LONG) 

DO 330 J=Jt1N, JNX 

IF ( <STt1L< J+2 )-STNL<L- 1 ) ) . LT . 8.*(STML(J+1 )-ST11L (L))) GO TO 320 
UR i TE < 6 , 3 1 0 ) XL CL). VL < L ) , ZL < L ) , XL C U+ 1 ) . YL C J+ 1 ) , ZL < J+ 1 ) 

3 1 0 FORMAT < 1 4H POINT DELETED . 1 0X . 3F 1 2 . 5 , 1 0X , 3F 12.5) 

STMLCL )=CRF*STMLCL >+STML< J+ 1 ) )/CRF+ 1 . ) 

XL CL )=CftF*XLCL )+XLC J+ 1 ) )/CAF+ 1 . > 

VLCL >=CP,F*VLCL )+VLCJ+ 1 ) >/CAF+ 1 . > 

ZL CL )= < AF+ZL C L )+ZL C J+ 1 ) ) / C AF+ 1 . > 

UXCL )=CAF*UXCL )+UXCU+ 1 ) )/CAF+ 1 . > 

UVCL )=CRF*UVCL )+UVCU+ 1 ) >/CRF+ 1 . ) 

IJZ CL )=CRF*U2CL )+IJZCU+ 1 ) >/CfiF+ 1 . ) 

GK 1 CL >=CAF*GK 1 CL )+GK 1 C J+ 1 > )/CAF+ 1 . ) 

GK2CL )=CAF*GK2CL )+GK2C J+ 1 ))/CAF+ 1 . ) 

H2CL )=CRF*H2CL )+H2CJ+ 1 ) )/CAF+ 1 . ) 

CP CL )= 1 .-(IJXCL >**2+lJVCL )**2+UZCL )**2 )/OSQ 
UflBSCL >=SQRT C 1 . -CPCL > ) 

AF=fiF+ 1 . 

GO TO 330 
320 RF= 1 . 

L=L+ 1 
K=J+ 1 

STMLCL )=STML GO 
XLCL)=XL(K) 

VLCL)=VLCK) 

ZLCL)=ZLCK) 

UXCL)=UXCK) 

UYCL )=UVCK) 
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UZ(L)=UZ<K) 

GKKL >=GK1(K> 

GK2CL >=GK2(K ) 

H2(L)=H2(K) 

CP(L)=CP(K) 

UflBS ( L )=UflBS < K ) 

NQUflDCL >=NQUflD(K ) 

330 CONTINUE 
L=L+ 1 

STML < L )=STML < JMflX > 

XL <L)=XL< JMflX) 

YL(L)=VL< JMflX) 

ZL(L)=ZL< JMflX) 

UXCL )=UX( JMflX ) 

UY<L)=UV( JMflX > 

UZCL)=UZ< JMflX) 

6KKL)=GKK JMflX) 

GK2<L>=GK2< JMflX) 

H2(L)=H2( JMflX) 

CP<L)=CP( JMflX) 

UflBS < L )=URBS < JMRX ) 

JMflX=L 

NQURDC JMflX >=NQlJflO< JMflX- 1 ) 

NQURDC JM I N )=NQUflD< JM ! N+ 1 ) 

1-iR ! TE <6,51 ) 

5 1 FORMAT < 4H0 I . 6X. , 1'riX . 9X , 1 HY . 

1 9X .1HZ, 09X , 2HUX , 8X , 2HUV , SX , 2HU2 , OOX , 

22HCP , SX , 2HK 1 , SX , 2HK2 . SX , 2HH2 , SX , 2HSL , SX , 1 HU , OX , 1 HP > 

IF CfiMflCH .ED. 0.) GOTO 1160 
C *+* COMPUTE COMPRESS I Bl LI TV CORRECTION 
DO 1150 J= JM I N , JMflX 

USD = < UX < J )**2+UV ( •- 1 >**2+UZ < J )+*2 ) /USQ 

USDS = USD 

SM = fiMflCH**2 

DO 1140 1=1,3 

F; = <1.+.2*SN*<1.-USDfl» 

IF CR .LT. .000001) R = .000001 
USDS = USD/R++2 . 5 
R - <1.+.2*SM*<1.-VSDB>> 

IF CR .LT. .000001) R = .000001 
USDC = US'D /F;**2 . 5 

1140 USDfi « CUSDC*USDfi-USDB**2 )/< VSDC+USDfl-2 . *USDB ) 

R = < I . + . 2*SM*( 1 . -USDfl ) ) 

IF CR .LT. .000001) R = .000001 
R = R** 1 . 25 
IJXCJ) = IJX< J )/F; 

UV<J) = UV< J)/R 
UZCJ) = UZCJ )/R 
UfiBSCJ) = SORT (USDfl) 

1 150 CPC J ) = <R**2 . 8- 1 . )/( . ?*SM > 

1150 CONTINUE 
K=0 

DO 53 I = JM I N, JMflX 
K=K+ 1 

53 HR I TE<6, 60 ) K, XL< I >, VL< I > . ZL< I ) , UX< I ) . UY< I >, UZC I ) . CP< I ) , 

1 GK 1C ! ) , GK2 < I ) , H2 ( ! ) , STML < I ) , UflBS Cl), NOUflD ( I 

60 FORMAT < IX, I3,3F10.5, 1X,3F10'.5. 1X.6F10.5, 16) 

8 FORMAT C3F 12.5) 

MR I TEC 17 > K, <XL< I ), YL< I ), ZL( I >, NOUflD < I >, I =JM I N, JMflX ) 
C I MR I TE .LE. 0 — WRITE SL,U,H2,K2 

C I WRITE . GE. 2 -- WRITE X,V,Z.CP 

C I NR i TE .EQ. 1 -- WRITE SL,U,H2,K2 AND X,Y,Z,CP 
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IF CiURiTE.GT. 1) GO TO 340 

MR I TE<7 > K, <STML( I >, UABSC I >, H2< I >, GK2< I ), I =JH I N , UMAX ) 
340 IF (IURITE.LT. 1) GO TO 400 

UR I TE<7 > K, <XL< I ), VLC I ), ZL< I ) , CPC I >, I =JM I N . JMRX > 

GO TO 400 

300 UR I TE (6 ,50) NSP(LL) 

WRITE <6, 65) 

GO TO 202 
400 CONTINUE 
00 TO 100 

REfiD NEXT SET OF STREAMLINES 

450 UP, ITEC6, 451 >IS,NQ 

451 FORMAT < 14H TAPE 04 ERROR, 2 1 4) 

550 ENDFILE 7 

REWIND 7 
ENDFILE 17 
REWIND 17 
REWIND 04 
STOP 7 
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- TRIRXIRL ELLlPSOiD 


INPUT 


FILE 




SAMPLE PROBLEM 
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280 5 150 150 


150 3 


.00001 0 
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0 
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0 0 


0 0 


1 . 00000 


. 00000 


.00000 


1 


1 


1 


0 


. ooooo 


.97861 


.00000 


. 10286 


1 


2 


1 


0 


. ooooo 


. 90789 


. 00000 


.20960 
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3 


1 


0 


. ooooo 


. 82360 


. 00000 


. 2S359 
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4 


1 


0 


.00000 


.71583 


.00000 


.34914 


1 


5 


1 


0 


. OOOOu 


. 62478 


. 00000 


. 39040 


1 


6 


1 


0 


. ooooo 


. 52537 


. 00000 


. 42544 


1 


“1 

f 


1 


0 


. ooooo 


.44721 


.00000 


.44721 


1 


8 


1 


0 


.00000 


. 37530 


. 00000 


. 45345 


1 


9 


1 


0 


. ooooo 


. 29822 


. 00000 


.47725 


1 


10 


1 


0 


. ooooo 


. 23692 


. ooooo 


. 48576 


1 


11 


1 


0 


. ooooo 


. 16967 


. ooooo 


. 49275 


1 


12 


1 

i 


o 

l_i 


. 00030 


. 1 1467 


. ooooo 


. 49670 


1 


13 


1 


0 


.00000 


. 05248 


. ooooo 


. 4993 1 


1 


14 


1 


0 


. ooooo 


. 00000 


. ooooo 


. 50000 


1 


15 


i 

l 


0 


. ooooo 


.99875 


. 10000 


. OOOOO 


2 


1 


1 


0 


. ooooo 


. 97739 


. 10000 


. 10273 


2 


2 


1 


0 


. ooooo 


. 90676 


. 10000 


. 20934 


2 


3 


1 


0 


. ooooo 


. 82257 


. 10000 


. 28323 


2 


4 


1 


o 


. 00000 


.71494 


. 10000 


. 34870 


2 


5 


1 


0 


. ooooo 


.62399 


. 10000 


. 3899 1 


2 


5 


1 


0 


. ooooo 


. 5247 1 


. 10000 


. 42490 


2 


7 


1 


0 


. ooooo 


.44655 


. 10000 


.44655 


2 


8 


1 


0 


. ooooo 


. 37483 


. 10000 


. 46287 


o 

j- 


9 


1 


0 


. ooooo 


. 29735 


. 10000 


. 47555 


2 


10 


1 


o 


ooooo 


. 23563 


. 10000 


.48516 


2 


11 


1 


o 


. ooooo 


. 16946 


. 10000 


.49213 


2 


12 


1 


0 


. ooooo 


. 1 1453 


. 10000 


. 49508 


z. 


13 


1 


0 


. ooooo 


.05241 


. 10000 


. 49859 


2 


14 


1 


0 


. ooooo 


. 00000 


. 10000 


. 4993? 


2 


15 


1 


0 


. ooooo 


99499 


. 20000 


. OOOOO 


2 


1 


1 


0 


. ooooo 


.97371 


.20000 


. 10234 


3 


2 


1 


0 


. ooooo 


. 90334 


. 20000 


. 20S55 


3 


3 


1 


0 


. ooooo 


.81947 


. 20000 


.28217 


3 


4 


1 


0 


.00000 


.71225 


.20000 


. 34739 


3 


5 


1 


0 


. ooooo 


.62164 


. 20000 


. 38845 




6 


1 


0 


.00300 


.52274 


. 20000 


.42330 


3 


7 


1 


0 


.00000 


. 44497 


. 20000 


.4449? 


3 


8 


1 


Ij 


. ooooo 


. 37342 


. 20000 


.45113 


3 


9 


1 


0 


. ooooo 


. 29672 


. 20000 


.47486 


•”« 


10 


1 


0 


. ooooo 


. 23574 


. 20000 


.48333 


o 


11 


1 


0 


. ooooo 


. 16882 


. 20000 


. 49028 




12 


1 


0 


. ooooo 


. 11410 


. 20000 


. 4942 1 


3 


13 


1 


0 


. ooooo 


. 05222 


. 20000 


.49581 


3 


14 


1 


0 


. ooooo 


.00000 


. 20000 


.49749 


O 

w 


15 


1 


0 


.00000 


. 98869 


.30000 


. OOOOO 


4 


1 


1 


0 


. ooooo 


. 95754 


. 30000 


. 10159 


4 


2 


1 


0 


. ooooo 


. 89762 


. 30000 


. 20723 


4 


3 


1 


0 


. ooooo 


.81428 


. 30000 


. 28038 


4 


4 


1 


0 


. ooooo 


. 70774 


.30000 


.34519 


4 


5 


1 


0 


. ooooo 


.61771 


. 30000 


. 38599 


4 


6 


1 


0 


. ooooo 


.51943 


.30000 


. 42052 


4 


( 


1 


0 


. ooooo 


.44215 


. 30000 


.44215 


4 


8 


1 


0 


. ooooo 


.37105 


. 30000 


. 4582 1 


4 


9 


1 


o 


. ooooo 


.29484 


.30000 


.47185 


4 


10 


1 


0 


. ooooo 


. 23424 


. 30000 


.4802? 


4 


1 1 


1 


0 


. ooooo 


. 16775 


.30000 


.48718 


4 


12 


1 


0 


. ooooo 



.000 .0 
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. 1 1338 
.05189 
. 00000 
. 97980 
. 95884 
. 88955 
.80696 
.70137 
,61215 
.51476 
.43818 
.36771 
.29219 
.23214 
. 16624 
. 1 1236 
.05142 
. 00000 
. 97980 
.95884 
. 88955 
. 80695 
.70137 
.61215 
.51476 
.43818 
. 3677 1 
.29219 
.23214 
. 15624 
. 1 1236 
.05142 
.00000 
. 96825 
.94754 
. 87905 
. 79745 
. 693 10 
.60494 
. 50869 
. 4330 1 
. 36338 
.28875 
.22940 
. 16428 
. 1 1 103 
0508 1 
.00000 
.95394 
.93354 
. 86507 
. 78566 
. 68286 
. 59600 
.50117 
. 42562 
. 3580 1 
.28448 
.22601 
. 16185 
. 10939 



. 30000 
.30000 
.30000 
. 40000 
.40000 
. 40000 
.40000 
.40000 
. 40000 
40000 
40000 
.40000 
40000 
. 40000 
.40000 
. 40000 
40000 
. 40000 
. 40000 
. 40000 
. 40000 
. 40000 
. 40000 
. 40000 
.40000 
. 40000 
. 40000 
.40000 
. 40000 
. 40000 
. 40000 
. 40000 
.40000 
. 50000 
. 50000 
. 50000 
.50000 
.50000 
. 50000 
. 5 C 030 
. 50000 
. 50000 
.50000 
. 50000 
. 5 C 000 
. 50000 
.50000 
. 50000 
. 60000 
. 60000 
. 60000 
. 60000 
. 60000 
. 60000 
.60000 
. 60000 
. 60000 
. 60000 
.60000 
. 60000 
. 60000 



.49108 
. 49366 
. 49434 
. 00000 
. 10078 
. 20537 
. 27786 
. 34208 
.38251 
.41684 
.43818 
. 45409 
.46761 
.47595 
. 48280 
. 48667 
. 48922 
.48990 
.00000 
. 10078 
. 20537 
. 27786 
. 34208 
.38251 
.41684 
.43818 
. 45409 
.46761 
. 47595 
.48280 
48667 
. 48922 
.48990 
. 00000 
. 09959 
. 20295 
. 27458 
. 33805 
.37801 
.41193 
.43301 
.44874 
. 46209 
. 47034 
.47710 
. 48093 
.48346 
48412 
. 00000 
.09812 
. 19995 
. 27052 
.33305 
.37242 
. 40584 
.42661 
.44211 
.45527 
. 46339 
. 47005 
. 47382 



4 13 

4 14 

4 15 

5 1 

5 2 

5 3 

5 4 

5 5 

5 6 

5 7 

5 8 

5 9 

5 10 

5 11 

5 12 

5 13 

5 14 

5 15 

6 1 
6 2 

5 3 

6 4 

6 5 

6 6 
6 7 

6 8 
6 9 

6 10 
6 1 1 
6 12 
6 13 

6 14 

6 15 

7 1 

7 2 

7 3 

7 4 

7 5 

7 6 

7 7 

7 8 

7 9 

7 10 

7 11 

7 12 

7 13 

7 14 

7 15 

8 1 

8 2 
y 3 
8 4 

8 5 

8 6 
8 7 

8 8 
8 9 

8 10 
8 11 
8 12 
8 13 



1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
1 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 C 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 



. 00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
00000 
. 00000 
.00000 
.00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
.00000 
.00000 
. 00000 
.00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 



.05005 
. 00000 
.93675 
91671 
. 85047 
.77150 
.67056 
.58525 
.49214 
.41893 
.35156 
. 27936 
.22194 
. 15894 
. 10742 
.04916 
.00000 
.91652 
. 8969 1 
.83210 
. 75484 
. 65507 
. 57262 
.48151 
. 40988 
. 34397 
. 27332 
.21714 
15550 
. 10510 
.04810 
. 00000 
.91552 
. 8969 1 
.83210 
. 75484 
. 65607 
.57262 
.48151 
. 40983 
. 34397 
. 27332 
.21714 
. 15550 
. 10510 
.04810 
. 00000 
.89303 
.87393 
.81077 
. 73550 
.63926 
.55794 
.45917 
. 39937 
.33515 
.25632 
.21158 
. 15152 
. 1024 1 
.04687 



.60000 
.60000 
.70000 
. 70000 
.70000 
. 70000 
. 70000 
.70000 
.70000 
. 70000 
. 70000 
.70000 
. 70000 
. 70000 
. 7G000 
. 70000 
. 70000 
. 80000 
.80000 
. 80000 
. 80000 
. 80000 
. 80000 
. 80000 
. 80000 
. 80000 
.80000 
. 80000 
. 80000 
. 8 ' JU ' J'J 
. 80000 
. 80000 
. 80000 
.80000 
. 80000 
.80000 
.80000 
.80000 
. 80000 
.80000 
.80000 
. 80000 
. 8D000 
.80000 
.80000 
. 80000 
. 80000 
. 90000 
. 90000 
. 90000 
. 90000 
. 90000 
. 90000 
.90000 
. 90000 
. 90000 
.90000 
. 90000 
. 90000 
.90000 
. 90000 



.47531 
. 47697 
.00000 
.09635 
. 19635 
.26565 
. 32705 
.36571 
.39853 
. 4 1893 
. 434 14 
.44706 
. 45504 
.46158 
. 46529 
. 46773 
.46837 
. 00000 
.09427 
. 19210 
.25991 
.31999 
. 3578 1 
. 38992 
.40988 
. 42476 
. 4374 1 
. 4452 1 
.45151 
4 55 23 
. 45763 
. 45826 
. 00000 
.09427 
. 19210 
.25991 
.31999 
.35781 
.38992 
. 40988 
.42476 
. 4374 1 
. 4452 1 
.45151 
45523 
. 45763 
. 45826 

. oocoo 

.09185 
. 18718 
.25325 
.31179 
. 34854 
.37993 
. 39937 
.41388 
. 42520 
. 43380 
. 44004 
,44357 
. 44590 



8 


14 


8 


15 


9 


1 


9 


2 


9 


3 


9 


4 


9 


5 


9 


6 


9 


7 


9 


8 


9 


9 


9 


10 


9 


11 


9 


12 


9 


13 


9 


14 


9 


15 


10 


1 


10 


2 


10 


■-« 

O 


10 


4 


10 


5 


10 


6 


10 


-> 

i 


10 


8 


10 


Q 


10 


10 


10 


11 


10 


12 


10 


13 


10 


14 


10 


15 


11 


1 


11 


2 


11 


3 


1 1 


4 


11 


5 


11 


6 


1 1 


7 


1 1 


8 


11 


9 


11 


10 


11 


11 


11 


12 


11 


13 


11 


14 


11 


15 


12 


1 


12 


2 


12 


3 


12 


4 


12 


5 


12 


6 


12 


7 


12 


8 


12 


9 


12 


10 


12 


11 


12 


12 


12 


13 


12 


14 



2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 

2 C 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 
2 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 

3 0 



. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
.00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. 00000 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 

. OOOOO 

. ooooo 
. ooooo 
. ooooo 
. ooooo 
.00000 
. ooooo 
. ooooo 

. OOOOO 
.00000 
. ooooo 
. oocoo 
. ooooo 
.00000 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
.00000 
.00000 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
.00000 
. ooooo 
. ooooo 
.00000 
. ooooo 
. ooooo 
. ooooo 
. ooooo 
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.00000 
.86603 
. 84750 
. 78626 
.71326 
.61993 
.54107 
. 45498 
.38730 
. 32502 
.25325 
.20518 
. 14594 
. 0993 1 
. 04545 

. ooooo 

.83516 
.81730 
. 75824 
. 68784 
.59784 
.52179 
. 43877 
. 37350 
.31343 
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